Guidelines for Remote Sensing of Evapotranspiration

This chapter describes the possibilities, the limitations and the future of remote sensing of evapotranspiration (ET). The principles behind the techniques of remote sensing of ET are presented systematically. The mathematical formulations of the key equations are used to highlight the critical parts and the variables that remote ET is most sensitive to. The focus will be on the input data. Which input data do we definitively need, and with what accuracy? How can we select the best methodology to estimate ET spatially? A number of new developments will be introduced, and priorities for the near future formulated. There is no global, validated ET product available today. We can find products of other components of the terrestrial water cycle, like rainfall and soil moisture, but not of ET. This means that remote estimation of ET is custom made, and that it requires specific skills. At first glance, this is surprising, because the idea of remote sensing of evapotranspiration is more than three decades old (Jackson et al., 1977; Jackson et al., 1987; Seguin, 1988). In this chapter we hope to clarify the reasons why the operational dissemination of remote sensing evapotranspiration products lags behind. The one fundamental problem with estimating ET is that it cannot be measured directly. This is well illustrated by borrowing an allegory from the evangelist Billy Graham: “I’ve seen the effects of the wind but I’ve never seen the wind”. This quote is certainly true, in a literal sense, for evapotranspiration. Evapotranspiration affects the water and energy balances, and it is these effects of evapotranspiration which are observed. For example, evapotranspiration reduces soil moisture content and its cools the land surface. By studying changes in soil moisture with lysimeters, or by studying patterns in land surface temperature with remote sensing, ET is estimated. The problem is that ET is not the only factor affecting the water and energy budgets of the surface. Other processes and physical properties play important roles as well. For this reason, quite a large number of input variables are needed to achieve reasonably accurate estimates of ET; estimates that are at least better than rule-of-thumb a priori estimates. This is the case for field techniques, and even more so for remote techniques. In remote sensing an additional issue is that the input data are not from just one source. Often data of different satellite platforms are needed in conjunction with meteorological data from ground stations. In the process of collecting the data, the user faces practical and scientific problems. The practical problems of the user include the data collection, merging the data in a GIS database, and programming the algorithm for the calculation of evapotranspiration. Each of


Introduction
This chapter describes the possibilities, the limitations and the future of remote sensing of evapotranspiration (ET).The principles behind the techniques of remote sensing of ET are presented systematically.The mathematical formulations of the key equations are used to highlight the critical parts and the variables that remote ET is most sensitive to.The focus will be on the input data.Which input data do we definitively need, and with what accuracy?How can we select the best methodology to estimate ET spatially?A number of new developments will be introduced, and priorities for the near future formulated.There is no global, validated ET product available today.We can find products of other components of the terrestrial water cycle, like rainfall and soil moisture, but not of ET.This means that remote estimation of ET is custom made, and that it requires specific skills.At first glance, this is surprising, because the idea of remote sensing of evapotranspiration is more than three decades old (Jackson et al., 1977;Jackson et al., 1987;Seguin, 1988).In this chapter we hope to clarify the reasons why the operational dissemination of remote sensing evapotranspiration products lags behind.The one fundamental problem with estimating ET is that it cannot be measured directly.This is well illustrated by borrowing an allegory from the evangelist Billy Graham: "I've seen the effects of the wind but I've never seen the wind".This quote is certainly true, in a literal sense, for evapotranspiration.Evapotranspiration affects the water and energy balances, and it is these effects of evapotranspiration which are observed.For example, evapotranspiration reduces soil moisture content and its cools the land surface.By studying changes in soil moisture with lysimeters, or by studying patterns in land surface temperature with remote sensing, ET is estimated.The problem is that ET is not the only factor affecting the water and energy budgets of the surface.Other processes and physical properties play important roles as well.For this reason, quite a large number of input variables are needed to achieve reasonably accurate estimates of ET; estimates that are at least better than rule-of-thumb a priori estimates.This is the case for field techniques, and even more so for remote techniques.In remote sensing an additional issue is that the input data are not from just one source.Often data of different satellite platforms are needed in conjunction with meteorological data from ground stations.In the process of collecting the data, the user faces practical and scientific problems.The practical problems of the user include the data collection, merging the data in a GIS database, and programming the algorithm for the calculation of evapotranspiration.Each of these steps requires expertise and access to dedicated GIS software packages, because most algorithms are not available 'on the shelf'.The scientific problem is that the procedure requires merging of data measured at different time and spatial scales.The merging of data for the calculation of evapotranspiration is a typical example of a data assimilation problem.Ideally, both the accuracy and the representativeness of each of the input data are taken into account in the calculation of the final product: the spatial map of evapotranspiration.Although algorithms for the calculation of evapotranspiration are available in the literature and on the internet, there is no consistent data assimilation procedure attached that calculates the accuracy and reliability of the final product.In recent years there have been a number of initiatives to build global products of evapotranspiration (Vinukullu et al., 2011;www.wacmos.org),addressing the above mentioned issues.The priorities for the near future are to establish a consistent way of merging the input data, improve the data assimilation techniques and to validate algorithms.In addition, new sources of data may be introduced.A promising tool is the use of satellite based laser altimetry for surface roughness estimation (Rosette et al., 2008).It is the aim of this chapter to focus on principles that the algorithms have in common, and on the input data.Reviews of the history of remote sensing of ET can be found in the scientific literature (Courault et al., 2005;Glenn et al., 2007;Gowda et al., 2007;Kalma et al., 2008).In addition to these reviews, we would like to provide some anchor points and guidelines for the selection of a methodology for estimating ET in basin hydrology.We will quantify and evaluate the error of each of the input data, and show how this error propagates into the final result.For this analysis we will use theoretical considerations, a remote sensing model, and a selection of field data.

Principles of remote sensing algorithms for evapotranspiration
Although remote sensing of evapotranspiration has evolved since the first initiatives in the 1970's, the fundamental principle has remained the same.All remote sensing based evapotranspiration estimates make use of the thermal and visible bands and the formulation of the energy balance of the surface.The instantaneous latent heat flux of evaporation is calculated as a residual of the energy balance, and this latent heat flux is in turn converted into an evapotranspiration rate after time integration.An inherent problem of this approach is that the errors in the various terms of the energy balance are affecting the latent heat flux in a manner that is difficult to predict.For this reason it is necessary to evaluate the different terms of the energy balance individually.In the evaluation of the remote sensing algorithms presented in this section, we will discuss the terms, and indicate at what spatial and temporal resolution the data can be collected.It will become clear that the land surface temperature is the most important state variable.It plays a crucial role in sensible heat flux, ground heat flux and the balance of long wave radiation.Apart from the collection of accurate land surface temperature data, important selection criteria for a methodology are the heterogeneity of the land cover, the topography and the spatial resolution (sampling) of the remote data.Neglecting the energy used in the process of photosynthesis, the instantaneous energy balance equation (EBE) over crops reads: R n is the net radiation remaining in the system, G the ground heat flux, H the sensible heat flux and E is the latent heat flux that is the energy consumed in evapotranspiration (all in W m -2 ).Radiation fluxes are positive when directed towards the land surface, the other fluxes are positive when pointed away from the surface.The partition of energy between the terms is largely controlled by the availability of water or moisture in the system.When moisture is not restricted, λE reaches a maximum and H is small.In order to estimate ET, Eq. 1 is solved for λE.When applied to remote sensor retrievals, R n is solved entirely from a combination of radiation counting at sensor level and few ground information.Ground or soil heat flux is a minor component in densely vegetated areas, but a large term in non-vegetated or sparsely vegetated areas (Heusinkveld et al., 2004).The importance of a better evaluation of the soil heat flux is gaining attention, mainly to ensure the EBE closure in such areas.The evaluation of H is the major difficulty.There are several models and approaches to solve for H (SEB models) and a number of parameters and assumptions are still under debate.The remote sensing models for ET mainly differ in the way H is treated.
In the following sections, the individual terms of the EBE (Eq. 1) will be discussed in further detail.A theoretical description is presented for each term in the EBE, in combination with a discussion on the feasibility of data acquisition from remote and ground sources.

Net radiation
Net radiation R n is the dominant term in the EBE, since it represents the source of energy that must be balanced by the thermodynamic equilibrium of the other terms.The net radiation can also be expressed as an electromagnetic balance of all incoming and outgoing radiation reaching and leaving a flat horizontal and homogeneous surface as: Where S is the shortwave radiation, nominally between 0.25 to 3m and L is the long wave radiation, nominally between 3 to 100 m.The arrows show the direction of the flux entering '' or leaving '' the system.Equation 2 is very convenient from the data acquisition point of view since each term can either be obtained from available models, or directly from instruments at ground stations or remote platforms.As remote sensors are positioned looking to Earth, they measure outgoing radiation only.The incoming fluxes must be either modelled or derived through alternative methodologies.
The instantaneous incoming shortwave radiation (also called global radiation), S, is commonly measured at ground stations by means of pyranometers or solarimeters.These instruments usually work in the shortwave broadband range (usually 0.305 -2.4 m).This range comprises almost 96% of the spectral interval of the solar irradiance.Recently there are remote sensing products and clearinghouses that account for the incoming and outgoing shortwave and long wave radiation.The use of them may reduce the need of permanently operational ground radiometers.
The outgoing shortwave radiation is the portion of the shortwave reflected back to the atmosphere.It is characterized by the albedo.The reflectance is the ratio between the reflected and the incoming radiation in a certain wavelength over an arbitrary horizontal plane.The integrated value over all visible bands defines the albedo, r 0 .Since albedo is a reflective property of the material, it can be evaluated from remote sensors multi-spectral bands, and the integration to full shortwave range is approached by a linear model that might include the atmospheric correction.The shortwave radiation balance reads: For all bodies, the total incident radiation is either reflected by the body, absorbed by it or transmitted through.This is expressed by the Kirchoff's law: where   is the reflectivity,   the transmissivity and   the absorptivity.A blackbody is defined as a body that absorbs all the radiation that receives.A blackbody is a physical abstraction that does not exist in nature.To keep a body temperature constant, it should emit the same radiation that absorbs.As a consequence a property of blackbodies, the absorptivity, is equal to the emissivity, and both are equal to 1, while reflectivity and transmissivity are equal to zero.Terrestrial materials behave more as grey bodies, meaning that part of the received radiation is reflected back to the atmosphere, or in other words, not all the energy that receives is absorbed.In order to keep the temperature constant, the absorbed radiation should equal the emission, so again emissivity is equal to absorptivity.Because the reflectivity is not zero, emissivity of real bodies is smaller than 1.
The longwave radiation terms are calculated with Planck's equation extended to real bodies.
A blackbody having a kinetic temperature T 0 [K] emits in a single wavelength a radiation that corresponds to: Where L  bb is the blackbody energy emission [W m -2 m -1 ] and  is the wavelength [m].
The kinetic temperature is the temperature as it would be measured by a standard thermometer in contact with the surface of the body.Emissivity ε λ at a chosen wavelength is the ratio of the radiation emitted by a real body at temperature T 0 to the radiation emitted by a blackbody at the same temperature.By definition, a blackbody has a constant emissivity equal to one for all wavelengths, whereas the real emissivity varies with wavelength.For natural bodies, the thermal emission can then be written as: Integration of L over all wavelengths leads to: where = 5.67 x10 -8 W m -2 K -4 is the Stefan-Boltzmann constant and  0 is a broadband surface emissivity.A remote sensor working within a spectral range of the thermal channels measures only a portion of L(T 0 ).The outgoing longwave radiation at any sensor channel is calculated by integration over the spectral range of the sensor: 00 () () The surface temperature T 0 is retrieved from Eq (8), once the surface emissivity   in the considered thermal channel is estimated.Once 'T 0 ' is obtained and  0 estimated, L is retrieved from Eq 7. Before the application of Eq (8), an atmospheric correction process is needed to derive L i sur (T 0 ) at the surface, because L i sat (T 0 ) as measured at the satellite sensor is affected by atmospheric interference.Atmospheric correction in the thermal range and in the shortwave is out of the scope of this chapter.We only mention that L i sur (T 0 ) can be obtained from L i sat (T 0 ) and using atmospheric correction model, in which water vapour and aerosol concentrations are the main input variables.
The incoming long wave radiation cannot be derived directly from remote sensors.It can either be determined from ground data or derived after atmospheric modelling.It varies with cloudiness (water vapour), air temperature and atmospheric constituents.For clear skies, the notion of effective thermal infrared emissivity of the atmosphere or apparent emissivity of the atmosphere ( ' a  ) introduces an overall emission value for all constituents.
If the air temperature T a at screen level is available, L is estimated as: There are several models simple to evaluate ' a  .The apparent emissivity of the atmosphere is usually estimated with equations based on vapour pressure and temperature at standard meteorological stations.For clear skies a common formulation, among others, is (Brutsaert, 1975): Where T a is the air temperature [K] and e a is the vapour pressure [mbar], everything measured at screen level.A portion L reaching the Earth surface is reflected back to the atmosphere.Since the surface is opaque the transmissivity is zero, the reflection of L↓ can be evaluated with Kirchoff law.As  0 describes the emissivity of a body in the thermal range, (1- 0 ) accounts for the reflection.The final expression for R n becomes: '4 '4 4 00 0 (1 Eq. 11 is valid for instantaneous observations.The conversion to a daily value is briefly discussed in Sect 2.4.It is not always necessary to carry out the calculations of Eqs 3-11 manually.Some organizations provide atmospherically corrected components of the radiation platforms directly.For example LandSaf (landsaf.meteo.pt)provides MeteoSat Second Generation (MSG) products of atmospherically corrected S and L with a 15 minute resolution and daily albedo for South America, Africa and Europe.An emissivity product will be released soon.Validation over ground based measurements for a site in Spain over sparse vegetation shows that these products are rather reliable (Fig 1 .).
As an alternative to the use of satellite data, a computation of the radiation terms from synoptic weather stations is also possible.The recommendations by the FAO (Doorenbos and Pruitt, 1977;Allen et al., 1998) could be followed.The daily short wave radiation S↓ day [MJm -2 day -1 ], is measured at agrometeorological stations with pyranometers and integrated to daytime hours.In most areas in the world, only sunshine hours are measured with periheliometers.In that case, the daily incoming shortwave radiation S↓ day can be obtained from the following empirical relationship:  The net daily shortwave radiation S day is estimated as in Eq. ( 3), assuming an average daily (sun hours only) albedo r 0day .The daily longwave radiation exchange between the surface and the atmosphere is very significant.Since on average the surface is warmer than the atmosphere and  0 > a ', there is usually a net loss of energy as thermal radiation from the ground.The daily net shortwave radiation L day [W m -2 ] between vegetation and soil on the one hand, and atmosphere and clouds on the other, can be represented by the following radiation law: '4 day ,da ,mean ( 273.15) Where ' a,day [-] is the daily net emissivity between the atmosphere and the ground, f a cloudiness factor and T a,mean is the mean daily air temperature at screen level [°C].Parameter ' a,day can be estimated from data from meteorological stations as: The cloudiness factor f is equal to 1 in case of a perfect clear day and 0 in a complete overcast day.In case the station has solar radiation data from pyranometers f can be calculated as:  (Maidment, 1992).
If only data on sunshine hours data are available, then:

Sensible heat flux
The sensible heat flux (H) is the exchange of heat through air as a result of a temperature gradient between the surface the atmosphere.Since the surface temperature during the day is usually higher than the air temperature, the sensible heat flux is normally directed upwards.During the night the situation may be reversed.Close to the surface, the sensible heat transport takes place mostly by diffusive processes, whereas at some distance away from the surface turbulent transport becomes more important.
The mathematical formulation of the sensible heat flux is based on the theory of mass transport of heat and momentum between the surface and the near-surface atmospheric environment (surface boundary layer).All existing remote sensing algorithms for turbulent sensible heat flux use the analogy of Ohm law of resistance driven by a gradient of temperature: where  a is the density of moist air [kg m -3 ], c p is the air specific heat at constant pressure [J kg -1 K -1 ], r ah is the aerodynamic resistance to heat transport between the surface and the reference level [s m -1 ] and T s -T a is the driving temperature gradient between the surface (with temperature T s ) and the reference height (with temperature T a ).Equation 18shows that the estimation of sensible heat flux has two main elements: a temperature difference between two heights and the corresponding resistance.As a first approximation we can conclude that the error in the sensible heat is linearly proportional to the error in the temperature gradient, and linearly proportional to the error in the inverse of the resistance.Equation 1shows that this error (W m -2 ) is directly transferred to the latent heat flux.We will now show that this is only approximately true, because the equation is not linear and the aerodynamic resistance itself depends on the temperature gradient.We will show that because of this, r ah can only be solved iteratively.
Understanding the physical concepts involved in the calculation of sensible heat flux, and in particular the aerodynamic resistance, is essential for an evaluation of remote sensing techniques.The evaluation of r ah is the most complicated issue of all in the whole EBE procedure for AET estimates.It is our experience that lack of or incomplete knowledge of the entire formulation, image pre-processing and atmospheric correction processes leads to severe flaws in the intermediate and final outputs.Many researchers are still seeking for alternatives, procedures and methods to improve the accuracy of ET estimates form the EBE -RS approach.The actual parameterization is not optimal in the sense that some sensitive information can only be strictly evaluated under controlled experimental research, and not in a routine fashion.
Near the ground two phenomena take place simultaneously in the transfer of heat between the surface and the atmosphere: free convection produced by temperature gradient T s -T a and forced convection by the dragging forces of the wind.Then, the estimation of the turbulent heat fluxes requires a description of the turbulent wind profile near the surface.The starting point of the analysis is the wind profile in a neutral atmosphere (no convection, and T s =T a ).In this situation and for an open site, the horizontal wind speed u [m s -1 ] varies logarithmically with height above the ground z [m]: B is usually replaced by A .ln (z om ) where z 0m is the aerodynamic roughness length of the surface for momentum transport and represents the value of z for which Eq 19 predicts In Eq 20, A must have the dimension of velocity and it should be independent of z since the profile description is given by the logarithmic term.The displacement height d is usually rather large: it ranges between 60 to 80 percent of the plant community height.Common values for well developed wheat are found in Verma and Bartfield (1979).A relation d= 0.67h is usually adopted for other vegetation types (Allen et al., 1998).An exact estimation can be carried out when a wind speed measurements at three or more heights are available.A plot of ln(z-d) versus wind speed should then give a straight line for the correctly calibrated value of d.Strictly speaking, d also depends on plant density; for sparse vegetation, d is often neglected.In many occasions, insufficient data are available to accurately predict its values for discrete crop canopies (Verma and Bartfield, 1979).If no valid field data are present, then it is suggested to leave d out of the equations altogether.The logarithmic wind profile in neutral conditions forms the basis for calculating the aerodynamic resistance for heat transport (and transport of water vapour) by mechanical turbulence.To include thermal turbulence (or convection or buoyancy), the Monin-Obukhov theory is needed in addition (Obukhov, 1946).We will now summarize the equations leading to the aerodynamic resistance of mechanical transport, followed by the modification for non-neutral conditions, when convection plays a role.The transfer of momentum in the direction of the flux takes place through molecular and turbulent eddy activity.Random vertical movements of the air cause air with different horizontal wind speeds to mix.This causes a momentum sink at the surface in a form of shear stress, .For convenience the shear stress is expressed as a function of a scalar u * , usually called eddy velocity or friction velocity [m s -1 ]: Since wind is produced by turbulent eddy motion, it is postulated that A in Eq 20 is proportional to the speed of the internal eddies.Then it can be demonstrated that: www.intechopen.comEvapotranspiration -Remote Sensing and Modeling 236 Where k is Von Kármán's constant, experimentally found to be 0.41.The final expression of the wind profile under neutral atmosphere is: Then the gradient of wind speed with height can be expressed as: In the theory of estimating the aerodynamic resistance from the wind profile, only vertical transport is considered.In the atmosphere the steepest gradients of heat, wind speed and humidity are found in the vertical direction.Horizontal variation is present in the order of tens of kilometres (Brutsaert, 2005), but these are considered negligible.This implies that horizontal advection effects (for example between pixels) are not considered in the remote sensing approach, a serious restrictions in patchy (wet and dry) environments.
In analogy to horizontal wind speed, heat and water vapour also have vertical profiles near the surface.Vertical mixing then causes a transport of heat and vapour too, resulting in vertical fluxes of sensible and latent heat.The three fluxes, of momentum (F u ), heat (F h ) and vapour (F v ), can be expressed as the covariance of vertical wind speed (w') and concentration of the admixture (u', T' and q'): These equations can be linked with the approach of electrical analogy (Eq 18) by approximating the covariance to simply the product of the vertical gradient of the quantities at two different heights (Brutsaert, 2005).A dimensionless parameter C is needed to fit the equality.This can be done for all three quantities.For example, for heat: In this case, C h will depend on the heights 1, 2, 3 and 4. It is convenient to choose the heights 4 and 3 equal to 2 and 1: Similarly, for momentum transport (shear stress): The coefficient C d can be calculated by combining Eqs 22, 24, 26 and 29.Using z 1 =z 0m , and considering that at z= z 0m , the wind speed is zero: In neutral conditions it can be assumed that C h =C v =C d , and thus: The appearance of the average temperature at height z 0m in Eq 31 is inconvenient.The roughness height, z 0h , changes with surface characteristics, atmospheric flow and thermal dynamic state of the surface (Blümel, 1999;Massman, 1999).It can be shown that: where B -1 is the inverse Stanton number, a dimensionless heat transfer coefficient.Free convection might alter the forced convective eddies generated by wind turbulence.
During daytime or when temperature decreases with height, convection amplifies the vertical eddy motions (unstable condition).During the night or when inversion conditions occur, and temperature increases with height, the horizontal eddy motions are enhanced (stable conditions).Mechanical turbulence and buoyancy coexists in a form of a hybrid regime known as mixconvection.Monin and Obukhov showed that these conditions eventually lead to an alteration of the wind and temperature profiles (Brutsaert, 1982).The Monin-Obukhov similarity theory uses dimensional analysis t o c o r r e c t t h e w i n d p r o f i l e p r o d u c e d b y buoyancy effects in such conditions.A non-dimensional correction factor for momentum transfer  m () is introduced to correct the wind profile gradient for conditions different from neutral, in which  is the ratio of thermal to mechanical turbulence: They introduced semi-empirical functions to correct the wind profile depending on the stability, based on dimensional analysis, of the form: www.intechopen.com

Evapotranspiration -Remote Sensing and Modeling 238
Where L is defined as the Monin-Obukhov length (L = z  ) [m], calculated as: Where g is the gravity constant (9.81 m s -2 ).Semi-empirical expressions for the stability corrections h  and m  can be found in the literature, for example Paulson (1970) and Brutsaert (1982).It is important to realize that L depends on air temperature and sensible heat flux, while sensible heat flux and air temperature in turn depend on L. For this reason, an iterative procedure is needed to calculate L, u * and H using Eqs 35-37.

Ground heat flux
The ground heat flux has received relatively little attention compared to the other terms.This is often justified, because ground heat flux is usually the smallest of all terms.Moreover, the 24-hour sum of ground heat flux is close to zero, because the heat absorbed during the day is released during the night.
At the moment of a satellite overpass, ground heat flux is not necessarily negligible.At midday it usually varies from 10% of net radiation for dense vegetation to 45% of net radiation for bare soil (Clothier et al., 1986).Often a vegetation cover dependent ratio between G and R n is assumed at satellite overpass (Kustas et al., 1990).
If more accurate estimates of ground heat flux are required, for example in areas with sparse vegetation, then remote estimates of ground heat flux are possible with the method of Van Wijk and De Vries (1963).For this method, diurnal cycles of land surface temperature and net radiation are needed (Verhoef, 2004;Murray and Verhoef, 2007); this means that time series of data of a geostationary satellite are required.An equation for ground heat flux can be derived the thermal diffusion equation, assuming a periodic land surface temperature: where  is the thermal inertia of the soil (J m -2 K -1 s -1/2 ), which depends on texture and soil moisture, t is time (s),  = (2/N) is the radial frequency (s -1 ), N the length of the time series [s], A and B integration coefficients [C], and n the number of harmonics.The coefficients A and B are fitted against the observed land surface temperature time series, for a chosen number of harmonics.The thermal inertia  can be estimated from soil texture and soil moisture, or calibrated against night time radiation, by assuming that night time radiation equals the night-time ground heat flux.

Latent heat flux
Latent heat flux is finally calculated as a residual of the energy balance (Eq.1).Because H, G and R n are instantaneous measurements, it is necessary to find a procedure to integrate to daily totals.A common way to carry out this integration, is by making use of the evaporative fraction, .The evaporative fraction (Brutsaert and Sugita, 1992) is the energy used for the evaporation process divided by the total amount of energy available for the evaporation process: It is assumed that the evaporative fraction remains constant throughout the day.

inst hrs
  (40) Assuming that the ground heat flux integrated over 24-hours is negligible, the evapotranspiration rate over 24 hours can be calculated as: Where = 2.0501-0.00236T water MJ kg -1 (T in C),  w = 1000 kg m -3 and R n24 is the average net The assumption of a constant evaporative fraction may lead to underestimates of daily evaporation, because the evaporative fraction in reality has a diurnal cycle with a concave shape (Gentine, et al., 2007).The concave shape is caused by changes in weather conditions (wind, advection, humidity), a phase difference between ground heat flux and net radiation, and stomatal regulation.There is an alternative to the assumption of constant evaporative fraction if hourly weather data are available.It may then be assumed that the ratio of actual to reference evaporation is constant over the day; hourly values of reference evaporation can be calculated (Allen et al., 2007).The ratio of actual to reference evaporation is more stable, because it eliminates the effects of diurnal variations in weather conditions.

Data requirements and sensitivity
Every remote sensing based SEB model requires a sequence of dedicated ground and remote sensing data to properly operate.Efforts increasingly focus on the remote estimation of the necessary variables, but ground data are still needed in addition.
All models require net radiation and land surface temperature retrieved from remote sensing.The additional required information varies among algorithms.As an example we list the input needed for the remote sensing model SEBS (Su, 2002).This model explicitly solves Eqs 35-37.It also includes an algorithm to estimate kB -1 from vegetation cover fraction.SEBS requires the following data, most of which cannot be retrieved from remote sensing, but is obtained from ground-based meteorological data instead: All meteorological input must be instantaneous information collected at the time of satellite overpass, interpolated and re-sampled to the pixel size.Other models require similar input.Two source models do not use the concept of kB -1 , but require separate resistances for soil and vegetation.Although the exact input data varies per algorithm, the most important are those related to the calculation of sensible heat flux, in particular the surface-air gradient and the corresponding aerodynamic resistance r ah .The success or failure of a SEB relies on the skills of the research team to extract realistic values for these two variables.For this reason, we will discuss these in more detail in the following sections.

The temperature gradient
For the temperature gradient we need to estimate both the air and the land surface temperature.The issue is that sensible heat flux is proportional to a difference between two temperatures which are obtained from two different sources in the same vertical.For this reason great care should be taken to retrieve both temperatures accurately.
For the air temperature at reference height, interpolated data of meteorological stations are commonly used.We need the air temperature well above the canopy, in the atmospheric surface layer, for which the aerodynamic resistance is defined.The standard measurement height in meteorological stations of 2 m cannot be used for vegetation taller than this height.Thus a conversion of temperatures from the meteorological stations to a higher reference height is needed.Another option is to use temperature profiles disseminated by organizations like EUMetsaf (www.eumetsat.int).
For the surface temperature it is necessary to take a closer look at the concepts first.As discussed before, the radiometric temperature is the temperature as it is retrieved from a remote radiometer by inverting Stefan-Boltzmann's law, assuming a bulk emissivity for the thermal spectrum range of the radiometer.The kinematic temperature is the real, contact temperature.A third definition is needed here: the aerodynamic temperature, which is hypothetic temperature obtained when extrapolating the vertical profile of air temperature to the depth z 0h .The aerodynamic temperature is a conceptual model parameter that is close to the kinetic temperature, but they are not equal.The reason is that kinematic temperature varies between the elements of the surface within a remote sensing pixel.For example, sunlit and shaded parts of the soil and canopy may have rather different temperatures.This is particularly the case in a heterogeneous landscape, where bare soil, vegetated and paved areas are mixed.It is even the case in a homogeneous land cover, where leaf temperatures may differ depending on their vertical position in the crown.This is illustrated in Fig 3 , showing the diurnal variations of contact temperatures of a needle forest in the Netherlands, the Speulderbos site, measured during a field campaign in The Netherlands on 16 June 2006 (Su et al., 2009).The lines are ensembles of 8 soil surface and 9 needle temperature sensors, mounted at different heights of the canopy of just a few trees.The heterogeneity of the soil and canopy temperatures will affect the radiometric surface temperature.The radiometric temperature is predominantly affected by the upper, visible, part of the canopy.Lower canopy layers also contribute to the outgoing upward radiation, but their contribution will be relatively low due to re-absorption of radiation.The radiometric temperature also depends on the solar angle and the observation angle of the satellite.We can see pronounced differences in the observed radiometric temperatures.Cleary visible is the hotspot, the situation where the solar zenith and azimuth angles equal those of the sensor.In the hotspot the radiometric temperature is higher than outside the hotspot.Radiometric temperatures are also higher at lower zenith angles compared to NADIR observations (vertically downward, in the centre of the graph).The differences in temperature are up to 2 C, indicating that care should be taken of the observation angle relative to the solar angle.It is also possible to exploit the differences in radiometric surface temperature observed at different angles in order to separate soil and canopy kinetic temperatures (Timmermans et al, 2009).How severe is an error of 2 C for the estimation of sensible heat flux?Equation 18shows that the error in sensible heat flux is proportional to the ratio of the temperature gradient to the resistance.This means that for the same error in the temperature gradient, the error in the sensible heat flux will be larger if the aerodynamic resistance is low than if the aerodynamic resistance is high.
In order to consider the sensitivity more precisely, we take the example of a situation where the aerodynamic resistance is low: the Speulderbos forest site in The Netherlands.This site is equipped with a 46-m tall eddy covariance measurement tower.Because of the low aerodynamic resistance, the sensitivity to temperature is expected to be relatively high.For this site, we calculated the friction velocity and the sensible heat flux with Eqs 35-37, using a canopy height of 30 m, the measured wind speed and temperature at 45 m height, radiometric temperature measured with a long-wave radiometer, and assuming that z 0m = 0.12 h and d = 0.67 h.accurate.During the night (stable conditions), the performance is worse, but this is not a large problem because the night-time sensible heat flux proved very small.However, there is a 50% error in afternoon sensible heat flux.Can this be related to an error in the surface temperature?Figure 6 shows the result of a sensitivity analysis to surface temperature.A consistent bias was added to the measured time series (x-axis), and the resulting root mean square error (RMSE) of friction velocity and sensible heat flux was calculated (y-axis).The RMSE reaches a minimum when surface temperature is 0.5 C above the measured value, but it rises to unacceptably high values of the absolute temperature bias is greater than 2 C.
In this example, field data of radiometric surface temperature were used.What if remote sensing data are available?Satellite products are available at either high temporal (geostationary satellites) or at high spatial resolution (polar orbiting satellites).Data are available at a spatial resolution of 3-5 km and a temporal resolution of 15 minutes (MeteoSat or GOES) to 1 km resolution at a daily time scale (AVHRR, MODIS or MERIS), or 60 m with a repetition time of weeks to months (LANDSAT, ASTER).The low temporal resolutions are not really useful, because of the dynamic nature of the turbulent heat fluxes.The daily revisits are useful provided that reasonable assumptions are made about the diurnal cycles of the fluxes (see Sect 2.4).The orbits are designed to overpass at the same solar time every day.The 15-minute intervals are ideal, but the spatial resolution makes the estimation of an effective aerodynamic resistance difficult, as we will see later.A final issue that needs to be considered is the topography.In areas with large elevation differences, the interpolation technique for air temperature data is crucial.Errors of several degrees in the air temperature are easily introduced if an incorrect adiabatic lapse rate is used.It is possible to circumvent the problem of estimating the temperature gradient by using an image based calibration (Bastiaanssen et al., 1998), in which assumptions are made for the energy balance state at the hottest and the coolest pixel in the image.In the first versions of this approach the calculated fluxes depended on the size of the image that was selected and the on assumption that the hottest pixel is dry, but more recent developments do not suffer from this drawback.The Mapping Evapotranspiration with Internalized Calibration (METRIC) model (Allen et al., 2007) uses reference evapotranspiration of alfalfa to calibrate the relation between the temperature gradient and the measured surface temperature.In the METRIC model it is assumed that the evaporation in the wettest pixel is 5% above the reference evapotranspiration, and the evaporation of the driest pixel is estimated with a soilvegetation-atmosphere model.This has the additional advantage that the evaporation values are bound to a realistic minimum and a realistic maximum rate.The METRIC model also accounts for topography by correcting radiation for slope and aspect and temperature for elevation using a local lapse rate.

Sensitivity to the aerodynamic resistance
The roughness length z 0m (and often displacement height is linked to it) is recognized as the main source of error in the remote estimate of ET.Currently, there are several methods that can be used to approach a good z 0m (see Table 2).When near surface wind speed and vegetation parameters (height and leaf area index) are available, the within-canopy turbulence model proposed by Massman (1999) can be used to estimate aerodynamic parameters, d, the displacement height, and, z 0m , the roughness height for momentum.This model has been shown by Su et al. (2001) to produce reliable estimates of the aerodynamic parameters.If only the height of the vegetation is available, the relationships proposed by Brutsaert (1982) can be used.If a detailed land use classification is available, for example based on LandSat images, the tabulated values of Wieringa (1993) can be used.By using literature values, errors in the canopy height of the order of decimetres to several metres are likely to occur, and errors in the roughness length in the order of decimetres.Sparse canopies require special attention.In sparse canopies the temperature differences between canopy and soil may be over 20 ºC.In addition, no canopy height can be defined, which makes it difficult to estimate the roughness length z 0m , and normally d is neglected.A solution to these problems is to program a two-source model (e.g.Norman et al., 1995).An alternative solution is to modify the parameter kB -1 to incorporate the differences in surface temperature implicitly in the value of z 0h (Verhoef et al., 1997).We will illustrate the latter solution with a simple example.The sparse canopy of our example is a study site in the province of León, Spain.The vegetation cover fraction is 11%, consisting of patches of 6-m tall Quercus ilex and Quercus pyrenaica.Data of an eddy covariance flux tower are used for validation of the satellite product.For this site, a roughness length of z 0m =0.2 m was assumed, and a displacement height of d=0.The friction velocity and sensible heat flux were again calculated from Eqs 35-37.For z 0h , a value of 0.02 m was initially assumed (kB -1 = 2.3), and for wind speed, the field measurements at the flux tower were used.For net radiation and surface temperature, 15minute interval MeteoSat Second Generation (MSG) satellite data were used.The top panels in Fig 8 show the results of the satellite based algorithm.The friction velocity observations are accurately reproduced, but the modelled sensible heat flux is extremely high, even double the net radiation.The overestimate is solved when we reduce z 0h by four orders of magnitude (kB -1 = 11.5).The reduction in z 0h needed to match the model with the observations is large.This problem was discussed earlier after the Hapex-Sahel measurement campaign (Verhoef et al., 1997).It was then concluded that the whole concept of kB -1 is questionable.It is indeed recommended to avoid the use of kB -1 , and this can be done in two ways: (1) by using more complicated two-source models for sparse vegetation, or (2) to use image-based calibration to relate surface temperature to a temperature gradient between two heights well above z 0h .The second approach is used in models such as SEBAL (Bastiaanssen et al., 1998) or METRIC (Allen et al., 2007).

Method
A model exists to estimate the kB -1 from vegetation density (Su et al., 2001).This model is used in the remote sensing algorithm SEBS (Su, 2002).However, care should be taken with any kB -1 model for areas where no detailed information on cover or other field data are available for calibration.
In the future, global maps of surface roughness may become timely available.Through synthesis of LiDAR with high resolution optical remote sensing, the roughness parameters have been successfully estimated spatially (Tian et al., 2011).Surface maps produced with laser satellites (NASA's ICESat and the future ICESat2) are also promising tools for estimating roughness (Roxette et al., 2008).

Conclusions
All remote sensing algorithms for ET make use of the energy balance equation (EBE).In this equation, latent heat flux is calculated as a residual of the energy balance.Net radiation can be estimated from remote sensing products relatively easily.Ground heat flux can only be retrieved with geostationary satellites for sparsely vegetated areas or bare land.It is usually a minor term in vegetated areas that causes relatively small errors in the final ET product.
The most critical component of the energy balance is the sensible heat flux.In the calculation of the sensible heat flux, both the temperature difference (land surface temperature minus the air temperature) and the aerodynamic resistance need careful attention.
In areas with high elevation differences, the errors in temperature are usually so high, and temperature correction using local lapse rates is necessary.In flat areas, a local sensitivity analysis is recommended.For forest, the accuracy of the temperature gradient should be better than 2 C in order to achieve reasonable results.In sparse vegetation two source models are preferred over single-source models, because in the latter, parameterization of z 0h on operational basis is no better than a wild guess.If a two-source model is not an option, then image based calibration using reference evaporation is a good alternative in these areas.Accurate roughness information (z 0m ) is required; the information is preferably verified and monitored on the ground.Satellite laser altimetry provides a promising tool for better roughness estimates in the near future.
a s is the fraction of the extraterrestrial radiation reaching the ground in a complete overcast day (when n=0), a s + b s the fraction of the extraterrestrial radiation reaching the ground in a complete clear day (n=N), n the duration of bright sunshine per day [hours], N the total daytime length [hours], S 0, day  is the terrestrial radiation [MJm -2 day -1 ].Local instrumentation can be used to calibrate a s and b s for local conditions.

Fig. 1 .
Fig. 1.Comparison between MeteoSat Second Generation radiation products (symbols) with 5-minute interval ground based measurements (lines) for a pixel with sparse vegetation in central Spain, for 5 July 2010.

Fig. 2 .
Fig. 2. Example plot of mean wind speed u against ln(z), measured above a 31 m tall forest in The Netherlands in August 2009.The intercept with the vertical axis leads to z 0m +d = 17.9 m.
1. Reference height z ref [m]: height from the ground where measurements of temperature, wind, pressure and specific humidity are made [m]. 2. Air Temperature at reference height (T a ) [°C]. 3. Specific humidity [kg.kg -1 ] or relative humidy [%], for calculation of emissivity of the sky.4. Wind speed at the reference height (u ref ) [m.s -1 ]. 5. Air Pressure at reference height [Pa].6.Air Pressure at land surface and reference height [Pa].7. The planetary boundary (PBL) height h i [m], required for the calculation of stability.It can be estimated by radiosounding or using atmospheric model outputs.By default h i =1000 m. 8.A map of vegetation heights, or alternatively, classes associated with vegetation height values, a map of Leaf Area Index (LAI) from which vegetation height is estimated.

Fig. 3 .
Fig. 3. Diurnal cycle of 8 soil and 9 needle contact temperature measurements (with NTC sensors) at the Speulderbos needle forest site in the Netherlands, on 16 June 2006.temperaturesof different elements in the canopy and the radiometric temperature.An example of an analysis carried out with SCOPE is shown in Fig 4.This figure shows simulated radiometric temperature of sparse but homogeneous crop as a function of the satellite observation azimuth (the counter-clockwise rotation angle from the top in the graph) and zenith angle (distance from the centre of the graph).We can see pronounced differences in the observed radiometric temperatures.Cleary visible is the hotspot, the situation where the solar zenith and azimuth angles equal those of the sensor.In the hotspot the radiometric temperature is higher than outside the hotspot.Radiometric temperatures are also higher at lower zenith angles compared to NADIR observations (vertically downward, in the centre of the graph).The differences in temperature are up to 2 C, indicating that care should be taken of the observation angle relative to the solar angle.It is also possible to exploit the differences in radiometric surface temperature observed at different angles in order to separate soil and canopy kinetic temperatures(Timmermans et  al, 2009).How severe is an error of 2 C for the estimation of sensible heat flux?Equation18shows that the error in sensible heat flux is proportional to the ratio of the temperature gradient to the resistance.This means that for the same error in the temperature gradient, the error in the sensible heat flux will be larger if the aerodynamic resistance is low than if the aerodynamic resistance is high.In order to consider the sensitivity more precisely, we take the example of a situation where the aerodynamic resistance is low: the Speulderbos forest site in The Netherlands.This site is equipped with a 46-m tall eddy covariance measurement tower.Because of the low aerodynamic resistance, the sensitivity to temperature is expected to be relatively high.For this site, we calculated the friction velocity and the sensible heat flux with Eqs 35-37, using a canopy height of 30 m, the measured wind speed and temperature at 45 m height, radiometric temperature measured with a long-wave radiometer, and assuming that z 0m = 0.12 h and d = 0.67 h.Figure5shows the results for 15-18 July 2009.The day-time friction velocity matches well with the measurements, showing that the calculation of aerodynamic resistance was

Figure 5 Fig. 5 .
Fig. 4. Hemispherical graph of simulated radiometric surface temperature of a thinned maize crop with a LAI of 0.25, as a function of viewing zenith angle and viewing azimuth angle (relative to the solar azimuth).Zenith angle varies with the radius, the azimuth angle (in italic) increases while rotating anticlockwise from north.The solar zenith angle was 48° (after Van der Tol et al., 2009).

Fig. 6 .
Fig. 6.Root mean square error of modelled sensible heat flux and friction velocity versus a forced bias in the observed radiometric surface temperatures for the Speulderbos forest site, for 14-19 July 2009.

Fig. 7 .
Fig. 7. Root mean square error of modelled sensible heat flux H and friction velocity u* velocity versus assumed canopy height of the Speulderbos forest site, for 14-19 July 2009.
Where a e is a correlation coefficient (ranging from 0.34 to 0.44, with a default of 0.34), b e a correlation coefficient (ranging from -0.14 to -0.25 with a default of -0.14), e d,mean the average vapour pressure at temperature[mbar].If true e d,mean is not available, then it can be calculated from daily average relative humidity RH mean and mean air temperature T a,mean [°C]:

Table 1 .
Typical values the coefficients a c , b c , a s and b s for arid and humid climates s , b s , a c and i c are calibration values to be estimated through specialized local studies which involve measuring longwave radiation values.Average values for a c and b c in arid and humid environments can be found in Table1:

Table 2 .
Su, 2002)or the estimation of z 0m (After:Su, 2002).When all of the above information is not available, then the aerodynamic parameters can be related to vegetation indices derived from satellite data.However in this case, care must be taken, because the vegetation indices saturate at higher vegetation densities and the relationships are vegetation type dependent.For example, characteristic of the land surface are sometimes calculated from indices like the Normalized Difference Vegetation Index (NDVI), but there is no reason why NDVI should have a universal relation with surface roughness.A grass field may have a similar NDVI to that of a forest, but a roughness length that is an order of magnitude smaller.For this reason, literature data or ground-truth data are indispensible for an accurate estimate of the surface roughness.A relation between NDVI and surface roughness can only be made for low vegetation, normally irrigated.In that case a non-linear relation with vegetation structure is first established by assigning a maximum value and a minimum value of height corresponding to values of NDVI.We illustrate the sensitivity of the aerodynamic resistance model with the data set of the Dutch forest site introduced in the previous section.Now, the height of the forest was varied between 5 and 45 m, and the RMSE of sensible heat flux and friction velocity evaluated (Fig 7).Note that the vertical scale in Fig 7 is much smaller than in Fig 6, which indicates that for forest, the model is less sensitive to errors in the canopy height than errors in surface temperature.