**A Distributed Benchmarking Framework for Actual ET Models**

## Yann Chemin

*International Water Management Institute Sri Lanka*

## **1. Introduction**

420 Evapotranspiration – Remote Sensing and Modeling

Gerosa G., Vitale M., Finco A., Manes F., Ballarin Denti A. and Cieslik S., 2005. Ozone

Grünhage L., Haenel H.D., Jager H.J., 2000. The exchange of ozone between vegetation and

Hicks B.B., Baldocchi, D.D., Meyers T.P., Hosker R.P., Matt D.R., 1987. A Preliminary

Holtslag A.A.M., van Ulden A.P., 1983. A simple scheme for daytime estimates of the

Jarvis P.G., 1976. The interpretation of the variations in leaf water potential and stomatal

Körner C., Scheel J., Bauer H., 1979. Maximum leaf diffusive conductance in vascular plants.

Leuning R., 1995. A critical appraisal of a combined stomatal photosynthesis model for C3

Marzuoli R., Gerosa G., Desotgiu R., Bussotti F., Ballarin-Denti A., 2008. Ozone fluxes and foliar

Monteith J.L., 1981. Evaporation and surface temperature. *Quarterly Journal of the Royal* 

Patwardhan S., Pavlick, R. Kleidon A., 2006. Does the empirical Ball-Berry law of stomatal

Schafer K.V.R., Oren R., Tenhunen J.D., 2000. The effect of tree height on crown level

Stewart J.B. (1988) Modelling surface conductance of pine forest. *Agricultural and Forest* 

Thom A.S., 1975. Momentum, mass and heat exchange of plant communities. In Vegetation

Unsworth M.H., Heagle A.S., Heck W.W., 1984. Gas Exchange in open field chambers – I.

Webb R.A., 1972. Use of the boundary line in the analysis of biological data. *Journal of* 

Measurement and analysis of atmospheric resistance to gas exchange. Atmospheric

stomatal conductance. *Plant, Cell & Environment* 23, 4 365-375.

and Atmosphere. Ed. J.L. Monteith. Academic Press, London.

*Populus berolinensis*): a dose–response analysis. Tree Physiology 29, 67-76. Monin A.S., Obukhov A.M., 1954. Basic laws of turbulent mixing in the atmosphere near the

environmental conditions. *Environmental Pollution* 152, 274–284.

*Environment* 39, 3255-3266.

22, 517–529.

*Environmental Pollution* 109, 373–392.

*Society of London*, Series B 273, 593–610.

plants. *Plant, Cell & Environment*. 18, 339–355.

TrudyGeofizich eskowo Instituta 151(24), 163–187.

*Union*, Fall Meeting 2006, abstract #H51C-0498.

Photosynthetica 13, 45-82.

*Meteorological Society* 107, 1–27.

*Meteorology* 43, 19–35.

*Environment* 18, 373–380.

*Horticultural Science* 47, 309-319.

quantities. *Water, Air and Soil Pollution* 36, 311–330.

and ozone exposure. An Open-Top Chambers experiment in South Alpine

uptake by an evergreen Mediterranean forest (*Quercus ilex*) in Italy. Part I: Micrometeorological flux measurements and flux partitioning. *Atmospheric* 

atmosphere: micrometeorological measurement techniques and models.

multiple resistance routine for deriving dry deposition velocities from measured

surface fluxes from routine weather data. *Journal of Climate and Applied Meteorology*

conductance found in canopies in the field. *Philosophical Transactions of the Royal* 

injury development in the ozone-sensitive poplar clone Oxford *(Populus maximowiczii x* 

ground. *Translation in Aerophysics of Air Pollution* (In: Fay, J.A., Hoult D.P. (Eds.), AIAA, New York, 1969, pp. 90–119). Akademija Nauk CCCP, Leningrad,

conductance emerge from maximization of productivity? *American Geophysical* 

With the various types of actual ET models being developed in the last 20 years, it becomes necessary to inter-compare methods. Most of already published ETa models comparisons address few number of models, and small to medium areas (Chemin et al., 2010; Gao & Long, 2008; García et al., 2007; Suleiman et al., 2008; Timmermans et al., 2007). With the large amount of remote sensing data covering the Earth, and the daily information available for the past ten years (i.e. Aqua/Terra-MODIS) for each pixel location, it becomes paramount to have a more complete comparison, in space and time.

To address this new experimental requirement, a distributed computing framework was designed, and created. The design architecture was built from original satellite datasets to various levels of processing until reaching the requirement of various ETa models input dataset. Each input product is computed once and reused in all ETa models requiring such input. This permits standardization of inputs as much as possible to zero-in variations of models to the models internals/specificities.

## **2. Theoretical points of observation**

## **2.1 Net radiation and soil heat flux**

In the two-source energy balance approach, like TSEB and SEBS differ from the single-source concept of SEBAL and METRIC in the sense that the radiation and energy balances have separate formulations for either bare soil or canopy. The energy balance at any instantaneous moment is expressed by equation Eq. 1:

$$Rn = G + H + LE\tag{1}$$

Where Rn is Net Radiation, G is soil heat flux, H is sensible heat flux and LE is latent heat of vaporization. This is what is appearing in single-source models like SEBAL and METRIC. Single source models concentrate on identifying Rn and G from astronomical and semi-empirical equations respectively, while H is being iteratively solved based on thermodynamically exceptional geographical locations, often referred in literature (Bastiaanssen, 1995) as *wet* and *dry* pixels, also the technique to identify them is referred in more recent literature as *end-members selection/identification* (Timmermans et al., 2007).

In two-source models, it is separated into bare soil and canopy energy balances as in Eq. 2 and 3, respectively:

$$\text{Rus} = \text{G} + \text{Hs} + \text{LEs} \tag{2}$$

Where Rns is the net radiation of bare soil surface, Hs is the sensible heat flux from bare soil, LEs is the latent heat of vaporization from soil surface.

$$R\mathfrak{nc} = H\mathfrak{c} + L\mathfrak{E}\mathfrak{c} \tag{3}$$

Where Rnc is the net radiation from canopy of crop, Hc is the sensible heat flux from canopy, LEc is the latent heat of vaporization of crop. Once the elements of those two equations are found, the fraction of vegetation cover (fc) is used to combine them into the area of a satellite remote sensing pixel, which is inherently a mixel of bare soil and canopy.

The Net Radiation is partitioned according to the formulation commonly used in two-sources model (Eq. 4 and 5), where the soil partition of Rn is an LAI-based extinction coefficient (Choudhury, 1989) with a coefficient C ranging from 0.3 to 0.7 (Friedl, 2002), depending on the arrangement of the canopy elements. Friedl (2002) mentions that a canopy with spherical (random) leaf angle distribution would lead to a C value of 0.5.

$$R\mathfrak{n}\mathfrak{s} = R\mathfrak{n}\,e^{\frac{-\zeta L A I}{\mathrm{Cov}(\mathrm{sumza})}}\tag{4}$$

$$Rn\mathcal{c} = Rn - Rns\tag{5}$$

Where LAI is the leaf area index, sunza is the sun zenith angle. Friedl (2002) mentions that he derived his soil heat flux formulation from his previous work (Friedl, 1996). It takes the already available soil fraction of net radiation and the cosine of the sun zenith angle (Eq. 6). A coefficient is then multiplied to those whereby soil type and moisture conditions are taken into consideration after (Choudhury et al., 1987).

$$\mathbf{G} = \mathbf{K} \mathbf{g} \, R \mathbf{n} \mathbf{s} \, \mathbf{C} \mathbf{s} (\mathbf{s} \,\mathrm{m} z \mathbf{a}) \tag{6}$$

Where Kg is the soil type and moisture condition coefficient in the soil heat flux. The Fraction of Vegetation cover is necessary to split the two-sources of heat transfer studied in such models. They are the soil surface (bare soil) and the vegetation canopy surface. The fraction of vegetation cover from Jia et al. (2003) quoting Baret et al. (1995) is developed as in Eq. 7:

$$f\mathcal{C} = 1 - [\frac{(NDVI - NDVI\_{min})}{(NDVI\_{min} - NDVI\_{max})}]^K\tag{7}$$

with K being taken as 0.4631 in Jia et al. (2003) and NDVImin at LAI=0 and NDVImax at LAI = +INF. As can be seen, a very large weight of potential deviation from the expected result is resting in the proper assessment of fc (Eq. 7). There are also uncertainties in the LAI raster input (Yang, Huang, Tan, Stroeve, Shabanov, Knyazikhin, Nemani & Myneni, 2006; Yang, Tan, Huang, Rautiainen, Shabanov, Wang, Privette, Huemmrich, Fensholt, Sandholt, Weiss, Ahl, Gower, Nemani, Knyazikhin & Myneni, 2006).

The soil heat flux computed for Bastiaanssen (1995), is what could be called a *partial contribution* of soil heat flux to the energy balance of the pixel, as the semi-empirical relationship is proportional to various elements of thermodynamic forcing within each pixel (Eq. 8).

$$G = \frac{Rn}{Albedo} \ T\_c \left( 0.0032 (\frac{Albedo}{r\_0}) + 0.0062 (\frac{Albedo}{r\_0})^2 \right) \left( 1 - 0.978 \text{NDVI}^4 \right) \tag{8}$$

with *Tc* the temperature in Celsius and *r*<sup>0</sup> the Albedo to apparent Albedo correction ranging 0.9 to 1.1 depending on the time of the day.

SEBS uses a two-source Albedo anchors stretching equation multiplied by the soil fraction of the pixel to extract a percentage of the net radiation as soil heat flux (Eq. 9).

$$\mathcal{G} = \text{Rn} \left( Albedo\_{dark} + (1 - f\_c) \left( Albedo\_{bright} - Albedo\_{dark} \right) \right) \tag{9}$$

Generic values are *Albedodark* = 0.05 and *Albedobright* = 0.35, while adjustements are made when concentrating on a specific land use, eventually.

#### **2.2 Monin-Obukhov Similarity Theory**

2 Will-be-set-by-IN-TECH

In two-source models, it is separated into bare soil and canopy energy balances as in Eq. 2

Where Rns is the net radiation of bare soil surface, Hs is the sensible heat flux from bare soil,

Where Rnc is the net radiation from canopy of crop, Hc is the sensible heat flux from canopy, LEc is the latent heat of vaporization of crop. Once the elements of those two equations are found, the fraction of vegetation cover (fc) is used to combine them into the area of a satellite

The Net Radiation is partitioned according to the formulation commonly used in two-sources model (Eq. 4 and 5), where the soil partition of Rn is an LAI-based extinction coefficient (Choudhury, 1989) with a coefficient C ranging from 0.3 to 0.7 (Friedl, 2002), depending on the arrangement of the canopy elements. Friedl (2002) mentions that a canopy with spherical

*Rns* <sup>=</sup> *Rn e* <sup>−</sup>*C LAI*

Where LAI is the leaf area index, sunza is the sun zenith angle. Friedl (2002) mentions that he derived his soil heat flux formulation from his previous work (Friedl, 1996). It takes the already available soil fraction of net radiation and the cosine of the sun zenith angle (Eq. 6). A coefficient is then multiplied to those whereby soil type and moisture conditions are taken

Where Kg is the soil type and moisture condition coefficient in the soil heat flux. The Fraction of Vegetation cover is necessary to split the two-sources of heat transfer studied in such models. They are the soil surface (bare soil) and the vegetation canopy surface. The fraction of vegetation cover from Jia et al. (2003) quoting Baret et al. (1995) is developed as in Eq. 7:

*f c* <sup>=</sup> <sup>1</sup> <sup>−</sup> [ (*NDV I* <sup>−</sup> *NDV Imin*)

with K being taken as 0.4631 in Jia et al. (2003) and NDVImin at LAI=0 and NDVImax at LAI = +INF. As can be seen, a very large weight of potential deviation from the expected result is resting in the proper assessment of fc (Eq. 7). There are also uncertainties in the LAI raster input (Yang, Huang, Tan, Stroeve, Shabanov, Knyazikhin, Nemani & Myneni, 2006; Yang, Tan, Huang, Rautiainen, Shabanov, Wang, Privette, Huemmrich, Fensholt, Sandholt, Weiss,

The soil heat flux computed for Bastiaanssen (1995), is what could be called a *partial contribution* of soil heat flux to the energy balance of the pixel, as the semi-empirical relationship is proportional to various elements of thermodynamic forcing within each pixel

) + 0.0062(

*Albedo r*0

*Albedo r*0

(*NDV Imin* − *NDV Imax*)

remote sensing pixel, which is inherently a mixel of bare soil and canopy.

(random) leaf angle distribution would lead to a C value of 0.5.

into consideration after (Choudhury et al., 1987).

Ahl, Gower, Nemani, Knyazikhin & Myneni, 2006).

*Albedo Tc* (0.0032(

*<sup>G</sup>* <sup>=</sup> *Rn*

(Eq. 8).

*Rns* = *G* + *Hs* + *LEs* (2)

*Rnc* = *Hc* + *LEc* (3)

*Cos*(*sunza*) (4)

*Rnc* = *Rn* − *Rns* (5)

*G* = *Kg Rns Cos*(*sunza*) (6)

]

*<sup>K</sup>* (7)

)2) (<sup>1</sup> <sup>−</sup> 0.978*NDV I*4) (8)

and 3, respectively:

LEs is the latent heat of vaporization from soil surface.

The Monin-Obukhov Similarity Theory (Monin & Obukhov, 1954) is being used in single source and two-source energy balance models. It is interesting to note that Monin & Obukhov (1954), in the development of their Monin-Obukhov Similarity Theory (MOST) considered the friction velocity to be about 5% of the geostrophic wind velocity having an average speed of 10m/s results in the friction velocity being around 0.5 m/s, and with the Coriolis parameter *l* = 10−4*s*−<sup>1</sup> and a tolerance of 20%, an estimate of the height of the surface layer is found at h=50m, that is also the DisALEXI blending height for air temperature (Norman et al., 2003). The dynamic velocity within this layer can be considered near to constant and the effect of Coriolis Force neglected (Monin & Obukhov, 1954). Under those conditions of neutral stratification the processes of turbulent mixing in the surface layer can be described by the logarithmic model of the boundary layer (Eq. 10).

$$\begin{aligned} L &= \frac{-1004\,\mu u^3 T}{kgH} \\ most\_x &= (1 - 16\,\frac{h}{L})^{\frac{1}{2}} \\ \psi\_{\text{fl}} &= 2\log(\frac{1 + most\_x^2}{2}) \\ \psi\_{\text{fl}} &= 2\log(\frac{1 + most\_x}{2}) + \log(\frac{1 + most\_x^2}{2}) - 2\,atan(most\_x) + 0.5\pi \end{aligned} \tag{10}$$

with *ψm*, *ψ<sup>h</sup>* the diabatic correction of momentum and heat through their changes of states, *mostx* a MOST internal parameter, *L* the Monin-Obukhov Length (MOL), *k* is the von Karman constant, g the gravity acceleration, *u* is the wind speed, *ρ* is the air density, *T* is the temperature and *h* is the height of interest (measurement height of the wind speed, roughness length, etc.).

Constraints to MOST as found in Bastiaanssen (1995) are of two types, first avoiding the latent heat flux input to be nil as its input location is in the denominator of the MOL equation (Equation Eq. 11), the second constraint is when the MOL is becoming positive, to force *ψ<sup>m</sup>* and *ψ<sup>h</sup>* to a ranged negative value (Bastiaanssen, 1995).

$$\begin{aligned} if(H=0.0): L &= -1000.0 \\ if(L>0.0): \psi\_h &= \psi\_m = -5\frac{2}{L} \end{aligned} \tag{11}$$

It turns out that Su (2002), extending the reach of his SEBS model to the GCM community has included a dual model for the convective processes within the Atmospheric Boundary Layer (ABL). Su (2002) followed the observations of Brutsaert (1999) that the ABL lower layer is either stable, either unstable and that the thickness of this lower layer is *α* = 10-15% of the ABL height, which is about *β* = 100-150 times the surface roughness. SEBS takes the highest from both as its estimation of *hst*, the height of ABL sublayer separation. If the reference height is lower than that, then the lower sublayer model is run, otherwise the upper sublayer models is used.

The cutline between the two sublayers of the ABL permits SEBS to process the lower layer (Atmospheric Surface Layer, ASL) under the MOST paradigm (Eq. 12), whether it proves unstable (generally in the day) or stable (generally in the night). The momentum and heat functions for the upper sublayer of the ABL where flow is laminar (free convection) can then be merged with the ASL by what Su (2002) calls a Bulk ABL Similarity (BAS) stability correction set of functions called here *ξ<sup>m</sup>* and *ξ<sup>h</sup>* for momentum and heat respectively (Eq. 13).

$$\begin{split} &if(\frac{-z\_0}{L}<0.0): \xi\_m = -2.2\operatorname{alog}(1+\frac{z\_0}{L})\\ &if(z\_0<\frac{a}{\beta}h\_i): \xi\_m = -\operatorname{alog}(a) + \psi\_m(\frac{-ah\_i}{L}) - \psi\_m(-\frac{z\_0}{L})\\ &if(z\_0\geq\frac{a}{\beta}h\_i): \xi\_m = \operatorname{alog}(\frac{h\_i}{\beta z\_0}) + \psi\_m(\frac{-\beta z\_0}{L}) - \psi\_m(-\frac{z\_0}{L}) \end{split} \tag{12}$$

$$\begin{split} &if(\frac{-z\_0}{L}<0.0): \xi\_h = -7.6\operatorname{alog}(1+\frac{z\_0}{L})\\ &if(z\_0<\frac{a}{\beta}h\_i): \xi\_h = -\operatorname{alog}(a) + \psi\_h(\frac{-ah\_i}{L}) - \psi\_h(-\frac{z\_0}{L})\\ &if(z\_0\geq\frac{a}{\beta}h\_i): \xi\_h = \operatorname{alog}(\frac{h\_i}{\beta z\_0}) + \psi\_h(\frac{-\beta z\_0}{L}) - \psi\_h(-\frac{z\_0}{L}) \end{split} \tag{13}$$

with *z*<sup>0</sup> the surface roughness for momentum, *hi* the height of ABL or Planetary Boundary Layer (PBL). The formulation of *ψ<sup>m</sup>* and *ψ<sup>h</sup>* functions are inherited from Beljaars & Holtslag (1991) and include either correction weights inside the standard equations in some cases (unstable conditions of *ψ<sup>m</sup>* and *ψh*), either a polynomial with exponential in other cases (stable conditions of *ψ<sup>m</sup>* and *ψh*). However, Beljaars & Holtslag (1991) stated categorically that the data described are characteristic for grassland and agricultural land *with sufficient water supply*.

#### **2.3 Roughness height**

Allen et al. (2005) mentions that METRIC and SEBAL do not require knowledge of crop type (no satellite based crop classification is needed). SEBAL relies on a type of semi-empirical equation relating NDVI to the roughness length (also called roughness height) for momentum and heat (Eq. 14).

$$z\_{0m} = e^{a + b \, NDVI} \tag{14}$$

Among many others, Chandrapala & Wimalasuriya (2003) proposed *a* = −5.5 and *b* = 5.8 for Sri Lanka using AVHRR NDVI images (sensor response curves, atmospheric correction and pixel size are all influencing NDVI response). Bastiaanssen (1995) preferred using extrema conditions to define *a* and *b*.

SEBS takes a ground truth added to mapping point of view and uses a look up table to translate land use raster maps into roughness length.

4 Will-be-set-by-IN-TECH

is either stable, either unstable and that the thickness of this lower layer is *α* = 10-15% of the ABL height, which is about *β* = 100-150 times the surface roughness. SEBS takes the highest from both as its estimation of *hst*, the height of ABL sublayer separation. If the reference height is lower than that, then the lower sublayer model is run, otherwise the upper sublayer models

The cutline between the two sublayers of the ABL permits SEBS to process the lower layer (Atmospheric Surface Layer, ASL) under the MOST paradigm (Eq. 12), whether it proves unstable (generally in the day) or stable (generally in the night). The momentum and heat functions for the upper sublayer of the ABL where flow is laminar (free convection) can then be merged with the ASL by what Su (2002) calls a Bulk ABL Similarity (BAS) stability correction set of functions called here *ξ<sup>m</sup>* and *ξ<sup>h</sup>* for momentum and heat respectively (Eq.

> *z*0 *L* )

> > −*αhi*

−*βz*<sup>0</sup>

*z*0 *L* )

−*αhi*

−*βz*<sup>0</sup>

*<sup>L</sup>* ) <sup>−</sup> *<sup>ψ</sup>m*(−*z*<sup>0</sup>

*<sup>L</sup>* ) <sup>−</sup> *<sup>ψ</sup>m*(−*z*<sup>0</sup>

*<sup>L</sup>* ) <sup>−</sup> *<sup>ψ</sup>h*(−*z*<sup>0</sup>

*<sup>L</sup>* ) <sup>−</sup> *<sup>ψ</sup>h*(−*z*<sup>0</sup>

*z*0*<sup>m</sup>* = *ea*+*b NDV I* (14)

*L* )

(12)

(13)

*L* )

*L* )

*L* )

*<sup>L</sup>* <sup>&</sup>lt; 0.0) : *<sup>ξ</sup><sup>m</sup>* <sup>=</sup> <sup>−</sup>2.2 *alog*(<sup>1</sup> <sup>+</sup>

*<sup>β</sup> hi*) : *<sup>ξ</sup><sup>m</sup>* <sup>=</sup> *alog*( *hi*

*<sup>L</sup>* <sup>&</sup>lt; 0.0) : *<sup>ξ</sup><sup>h</sup>* <sup>=</sup> <sup>−</sup>7.6 *alog*(<sup>1</sup> <sup>+</sup>

*<sup>β</sup> hi*) : *<sup>ξ</sup><sup>h</sup>* <sup>=</sup> *alog*( *hi*

*<sup>β</sup> hi*) : *<sup>ξ</sup><sup>h</sup>* <sup>=</sup> <sup>−</sup>*alog*(*α*) + *<sup>ψ</sup>h*(

*<sup>β</sup> hi*) : *<sup>ξ</sup><sup>m</sup>* <sup>=</sup> <sup>−</sup>*alog*(*α*) + *<sup>ψ</sup>m*(

*βz*<sup>0</sup>

*βz*<sup>0</sup>

with *z*<sup>0</sup> the surface roughness for momentum, *hi* the height of ABL or Planetary Boundary Layer (PBL). The formulation of *ψ<sup>m</sup>* and *ψ<sup>h</sup>* functions are inherited from Beljaars & Holtslag (1991) and include either correction weights inside the standard equations in some cases (unstable conditions of *ψ<sup>m</sup>* and *ψh*), either a polynomial with exponential in other cases (stable conditions of *ψ<sup>m</sup>* and *ψh*). However, Beljaars & Holtslag (1991) stated categorically that the data described are characteristic for grassland and agricultural land *with sufficient water supply*.

Allen et al. (2005) mentions that METRIC and SEBAL do not require knowledge of crop type (no satellite based crop classification is needed). SEBAL relies on a type of semi-empirical equation relating NDVI to the roughness length (also called roughness height) for momentum

Among many others, Chandrapala & Wimalasuriya (2003) proposed *a* = −5.5 and *b* = 5.8 for Sri Lanka using AVHRR NDVI images (sensor response curves, atmospheric correction and pixel size are all influencing NDVI response). Bastiaanssen (1995) preferred using extrema

SEBS takes a ground truth added to mapping point of view and uses a look up table to

) + *ψm*(

) + *ψh*(

is used.

13).

*i f*( −*z*<sup>0</sup>

*i f*(*z*<sup>0</sup> <

*i f*(*z*<sup>0</sup> ≥

*i f*(*z*<sup>0</sup> <

*i f*(*z*<sup>0</sup> ≥

**2.3 Roughness height**

conditions to define *a* and *b*.

and heat (Eq. 14).

*i f*( −*z*<sup>0</sup> *α*

*α*

*α*

*α*

translate land use raster maps into roughness length.


Table 1. Boundary conditions on *z*0*<sup>m</sup>* from NDVI satellite data after Bastiaanssen (1995)

While SEBAL and METRIC use a fixed conversion rate between roughness length for momentum and roughness length for heat, Su (2002) is introducing in SEBS the use of the exponential of *kB*−<sup>1</sup> (Eq. 15).

$$\begin{aligned} z\_{0m} &= 0.136 \, h\_{\text{vogetation}}\\ z\_{0h} &= \frac{z\_{0m}}{e^{kB^{-1}}} \end{aligned} \tag{15}$$

with *hvegetation* the height of the vegetation relating to the roughness length for momentum *z*0*<sup>m</sup>* from Brutsaert (1982), *z*0*<sup>h</sup>* the roughness length for heat, Su (2002) refers to the work of Massman (1999) on the combined *von Karman constant - sublayer Stanton number (kB*−1*)*, where *B*−<sup>1</sup> is defined by Gieske (2007) as in Eq. 16:

$$\begin{aligned} B^{-1} &= St\_k^{-1} - \mathbb{C}\_d^{-\frac{1}{2}} \\ \frac{St\_k^{-1}}{\mu\_\*} &= \frac{\rho \mathbb{C}\_p \Delta T}{H} \end{aligned} \tag{16}$$

*Cd* is the drag coefficient, *St*−<sup>1</sup> *<sup>k</sup>* is the roughness Stanton number, *u*<sup>∗</sup> is the friction velocity, *ρ* the air density, *Cp* the specific heat and ΔT the temperature difference.

#### **2.4 Aerodynamic roughness for heat**

The aerodynamic resistance (roughness) for heat *rah* is an input to the sensible heat flux and is often a source of concentration in ET models based on energy balance. For logical reasons, as the parameterization of *rah* needs prior knowledge of the state of the sensible heat flux to enable knowledge of the MOL to parameterize *ψ<sup>m</sup>* and *ψ<sup>h</sup>* the diabatic correction of momentum and heat through their changes of states. In turn, *ψ<sup>m</sup>* and *ψh*, offset the logarithmic relation of the observation height to the respective roughness lengths (*z*0*<sup>m</sup>* and *z*0*h*) being the driving force to curve the relation to the wind shear profile.

SEBS resistance at the wet limiting case *rewet* (Eq. 17) is using the MOL as *Lwet* configured to use all the energy available for evaporation, which in turn is used in conjunction with the reference height (*zref* ) either in the computation of the *ψ<sup>h</sup>* or *χ<sup>h</sup>* whether ASL or BAS models are at work.

$$re\_{\text{net}} = \frac{alog(\frac{z\_{\text{ref}}}{L\_{\text{net}}} - [\psi|\chi]\_{\hbar})}{ku^\*} \tag{17}$$

#### **2.5 Ground to air temperature difference**

Allen et al. (2005) mentions that METRIC is a variant of the important model SEBAL and that it has been extended to provide tighter integration with ground-based reference ET. SEBAL formulation for the ground to air temperature difference (*dT*) is estimated (Eq. 18) as an affine function with two extreme conditions found in the satellite image processed (Bastiaanssen, 1995).

$$\begin{aligned} pixel\_{coll} &\mapsto dT = 0\\ pixel\_{hot} &\mapsto dT = \frac{(Rn - G) \times r\_{all}}{\rho \mathbb{C}\_p} \\ dT &= a + b \times T \end{aligned} \tag{18}$$

with *rah* the aerodynamic resistance for heat, *pixelcold* and *pixelhot* the end-members defined in Bastiaanssen (1995); Bastiaanssen et al. (1998); Timmermans et al. (2007) and that are the *signet ring* of the model.

METRIC formulation (Eq. 19) includes the reference ET paradigm found in Allen et al. (1998), the *Kc* × *ETo* crop ET, also called ETc, being translated into METRIC as *k* × *ETr* (Allen et al., 2005), practically using Alfalfa at full growth as anchor point in their Idaho study. Some extra metorological data being available when using METRIC, permits a daily surface soil water balance to be run to enforce conditions on the *dry/hot* pixel energy balance and effective dT. Selection of the extreme pixels is focused on cropped area as much as possible.

$$\begin{aligned} pixel\_{cold} &\mapsto dT = \frac{(Rn - G - k \times ET\_r) \times r\_{ah}}{\rho \mathbb{C}\_p} \\ pixel\_{hot} &\mapsto dT = \frac{(Rn - G) \times r\_{ah}}{\rho \mathbb{C}\_p} \\ dT &= a + b \times T \end{aligned} \tag{19}$$

SEBS (Su, 2002) is computing dT (Eq. 20) from a surface skin virtual temperature (*T*0) and PBL virtual temperature (*Tpbl*).

$$\begin{aligned} T\_v &= \frac{\log(h\_{pbl} - h\_{disp})}{\log(h\_u - h\_{disp})} \\ T\_{pbl} &= \frac{T\_s \times (1 - f\_c) + T\_v \times f\_c}{\frac{1 - DEM}{44331.0}^{1.5029}} \\ T\_0 &= \frac{T\_c}{\frac{1 - DEM}{44331.0}^{1.5029}} \end{aligned} \tag{20}$$
 
$$dT = T\_0 - T\_{pbl}$$

Where *Ts* is the soil temperature, *Tv* is the (virtual) vegetation canopy temperature, *fc* is the fraction of vegetation cover, *DEM* is the elevation and *Tc* is the satellite sensed temperature in Celsius, *hu* is the wind speed measurement height. *hdisp* is the displacement height being 0.65 of the canopy height in SEBS, Monin & Obukhov (1954) mention that for observations made at height superior to 1 meter, the displacement height can be nullified. The blending height (*hpbl*) is given by an external mean or if no data is available, default value of 200m is used (same as in SEBAL) or 1000m.

#### **2.6 Actual ET**

Single sources models (METRIC and SEBAL) have promoted a particular way of closing the energy-balance (Eq. 21), using a ratio of fluxes called the evaporative fraction (Λ; Eq. 21).

6 Will-be-set-by-IN-TECH

*pixelhot* �→ *dT* <sup>=</sup> (*Rn* <sup>−</sup> *<sup>G</sup>*) <sup>×</sup> *rah*

with *rah* the aerodynamic resistance for heat, *pixelcold* and *pixelhot* the end-members defined in Bastiaanssen (1995); Bastiaanssen et al. (1998); Timmermans et al. (2007) and that are the

METRIC formulation (Eq. 19) includes the reference ET paradigm found in Allen et al. (1998), the *Kc* × *ETo* crop ET, also called ETc, being translated into METRIC as *k* × *ETr* (Allen et al., 2005), practically using Alfalfa at full growth as anchor point in their Idaho study. Some extra metorological data being available when using METRIC, permits a daily surface soil water balance to be run to enforce conditions on the *dry/hot* pixel energy balance and effective dT.

*pixelcold* �→ *dT* <sup>=</sup> (*Rn* <sup>−</sup> *<sup>G</sup>* <sup>−</sup> *<sup>k</sup>* <sup>×</sup> *ETr*) <sup>×</sup> *rah*

SEBS (Su, 2002) is computing dT (Eq. 20) from a surface skin virtual temperature (*T*0) and PBL

*Tpbl* <sup>=</sup> *Ts* <sup>×</sup> (<sup>1</sup> <sup>−</sup> *fc*) + *Tv* <sup>×</sup> *fc* 1−*DEM* 44331.0

1.5029

Where *Ts* is the soil temperature, *Tv* is the (virtual) vegetation canopy temperature, *fc* is the fraction of vegetation cover, *DEM* is the elevation and *Tc* is the satellite sensed temperature in Celsius, *hu* is the wind speed measurement height. *hdisp* is the displacement height being 0.65 of the canopy height in SEBS, Monin & Obukhov (1954) mention that for observations made at height superior to 1 meter, the displacement height can be nullified. The blending height (*hpbl*) is given by an external mean or if no data is available, default value of 200m is used

Single sources models (METRIC and SEBAL) have promoted a particular way of closing the energy-balance (Eq. 21), using a ratio of fluxes called the evaporative fraction (Λ; Eq. 21).

1.5029

*ρCp*

*ρCp*

*ρCp*

(18)

(19)

(20)

*pixelcold* �→ *dT* = 0

*dT* = *a* + *b* × *T*

Selection of the extreme pixels is focused on cropped area as much as possible.

*pixelhot* �→ *dT* <sup>=</sup> (*Rn* <sup>−</sup> *<sup>G</sup>*) <sup>×</sup> *rah*

*Tv* <sup>=</sup> *log*(*hpbl* <sup>−</sup> *hdisp*) *log*(*hu* − *hdisp*)

*<sup>T</sup>*<sup>0</sup> <sup>=</sup> *Tc* 1−*DEM* 44331.0

*dT* = *T*<sup>0</sup> − *Tpbl*

*dT* = *a* + *b* × *T*

1995).

*signet ring* of the model.

virtual temperature (*Tpbl*).

(same as in SEBAL) or 1000m.

**2.6 Actual ET**

$$
\Lambda = \frac{Rn - G - H}{Rn - G} \tag{21}
$$

This instantaneous ratio (could be called an efficiency, since it is unitless) is multiplied with a potential value of ET in order to provide with an Actual ET (equation Eq. 22).

$$ET\_a = \Lambda \, ET\_{potential} \tag{22}$$

METRIC (Allen et al., 2007) evaporative fraction is used with the potential ET based on Penman-Monteith method as published in Allen et al. (1998) in order to produce an actual ET estimation. In contrast SEBAL (Bastiaanssen, 1995; Bastiaanssen et al., 1998) evaporative fraction is used with the potential ET computed from exo-atmospheric solar radiation, a single-way atmospheric transmissivity, and the reflected proportion from Albedo.

SEBS (Su, 2002) is computing an instantaneous latent heat flux (*LE*) from the evaporative fraction and the subtraction of net radiation and soil heat flux (Eq. 23). However, using what is called in SEBS the *relative evaporation*, the pixel value of *H* the sensible heat flux is allowed to vary only between *Hwet* and *Hdry*.

$$\begin{aligned} H\_{dry} &= Rn - G\\ H\_{wet} &= \frac{(Rn - G) - \frac{(\frac{\mu C\_p}{r\_{\text{at}}} \varepsilon\_{\text{at}})}{\gamma}}{\frac{1 + slope}{\gamma}}\\ E\_{wet} &\leq H\_{wet} \leq H \leq H\_{dry} \\ E\_{relative} &= 1 - \frac{H - H\_{wet}}{H\_{dry} - H\_{wet}}\\ \Lambda &= E\_{relative} \frac{Rn - G - H\_{wet}}{Rn - G} \end{aligned} \tag{23}$$

*Hwet* is derived from the Priestley-Taylor equation (Priestley & Taylor, 1972), from which *γ* and *slope* come from, *re*\_*wet* is the resistance at the *wet* limiting case. *Erelative* is the relative evaporation.

SSEB (Senay et al., 2007) is evaluating a regional approximation of the evaporative fraction Λ as a ranging of satellite-based temperature products (Eq. 24) and uses the reference ET found in Allen et al. (1998) as a mean to compute the actual ET.

$$
\Lambda = \frac{T\_{hot} - T}{T\_{hot} - T\_{cold}} \tag{24}
$$

#### **3. Methodology**

#### **3.1 Conceptual design and processing flow**

The distributed framework is a Linux system based on GDAL library (GDAL, 2011) and C programming, enhanced with a distributed language called OpenMP (OpenMP, 2011), used essentially for data distribution as seen in (Chemin, 2010).

As the conceptual architecture of the framework (Fig. 1) signifies, there are several layers of processing involved. Initially, downloaded satellite imagery is located in a single directory (referred as *RS data* in Fig. 1).

During the pre-processing phase, a parsing agent will select images of the same kind and same date and stitch them together, then reproject them to the projection system selected. Upon completion, a second agent will perform a relatively conservative quality check and assign null value to failed pixels according to the satellite imagery information available.

Fig. 1. Architecture Concept of the Framework

Once these two steps are performed, the raster maps are tagged as *products*. Those products will be shared in between any of the ET models that require such type of input.

Some ET models require some higher-level input raster maps, by this, we define *higher-level product* as a raster that requires at least one *product* as defined above as precursor to its creation.

#### **3.2 Meteorological data**

Meteorological data (referred as *Point data* in Fig. 1) is encoded with Fourier Transforms (FT) in function of cumulative day of year from the beginning of the satellite imagery data set. This insures both faithfulness of the data and high-portability as well as an elegant way to summarize a complex and variable non-spatial dataset.

Actual state of the research in meteorological data time-series encoding for this framework is to develop an array of geo-tagged FT. Also under consideration are Wavelet Transforms (WT), for reasons of time tagging. This array will be used to interpolate on-the-fly from an appropriate number of neighbours and with an interpolation algorithm suiting best the operator requirement. The reason for this, is that it transfers the load from storage requirements (which can be heavy for daily meteorological raster datasets) to thread computing that is benefiting from distributed speed-up as only a marginal number of equations will be added by such process compared to the ET models processing load.

## **4. Implementation**

8 Will-be-set-by-IN-TECH

During the pre-processing phase, a parsing agent will select images of the same kind and same date and stitch them together, then reproject them to the projection system selected. Upon completion, a second agent will perform a relatively conservative quality check and assign null value to failed pixels according to the satellite imagery information available.

Once these two steps are performed, the raster maps are tagged as *products*. Those products

Some ET models require some higher-level input raster maps, by this, we define *higher-level product* as a raster that requires at least one *product* as defined above as precursor to its creation.

Meteorological data (referred as *Point data* in Fig. 1) is encoded with Fourier Transforms (FT) in function of cumulative day of year from the beginning of the satellite imagery data set. This insures both faithfulness of the data and high-portability as well as an elegant way to

Actual state of the research in meteorological data time-series encoding for this framework is to develop an array of geo-tagged FT. Also under consideration are Wavelet Transforms (WT), for reasons of time tagging. This array will be used to interpolate on-the-fly from an appropriate number of neighbours and with an interpolation algorithm suiting best the operator requirement. The reason for this, is that it transfers the load from storage requirements (which can be heavy for daily meteorological raster datasets) to thread computing that is benefiting from distributed speed-up as only a marginal number of

equations will be added by such process compared to the ET models processing load.

will be shared in between any of the ET models that require such type of input.

Fig. 1. Architecture Concept of the Framework

summarize a complex and variable non-spatial dataset.

**3.2 Meteorological data**

Practically (Fig. 2), MODIS datasets are grouped by products and by day and batch processed each in one core of the computer in parallel. This involves format changing, merging tiles, reprojecting, renaming outputs according to the nomenclature of the processing system. The tools involved in that step are either standard Linux Shell tools, either part of GDAL (2011) standard tools (i.e. *gdalwarp* and *gdal\_translate*). Both of these tools are still essentially sequential programs at this time, thus, they are being sent to each core in a distributed manner through the Shell with a check loop to ensure that there is at all time the same number of programs running as there are cores/threads available in the CPU architecture. It becomes clear that for each new leap in number of cores in future commercial offerings, the framework will automatically increase its processing capacity to the new enlarged number of cores/threads available, thus also reducing by the same factor the time needed to process a given number of satellite images.

Fig. 2. Architecture Implementation of the Framework

Models that are already inside the framework account to SSEB from Senay et al. (2007), METRIC from Allen et al. (2007), SEBAL from Bastiaanssen et al. (1998) using the work from Alexandridis et al. (2009), in progress are SEBS from Su (2002) and TSEB from both Kustas & Norman (1999) and Norman et al. (1995).

One of the many candidate for inclusion was the Two-Source Algorithm (TSA) from Yunhao et al. (2005). After extensive calculus and referring to the mathematic academia, it became clear that both temperature equations below when combined to extract *Tvegetation* and *Tsoil* (Eq. 25) have a large amount of solutions even within the constraining dimension of surface skin temperature range. The available information to possibly close the equations unknown parameters are relating to satellite *T* after Price (1984) and the satellite-based emissivity (*�*), partitioned into *�vegetation* = 0.93 & *�soil* = 0.97 quoting Sui et al. (1997). Unless missing information or other sources of solution dimension constraints are published, the proposition of the model is so far considered impossible and not included in the framework until such condition is being met through publication.

$$\begin{aligned} T &= f\_{\text{c}} \, T\_{\text{vegetation}} + (1 - f\_{\text{c}}) \, T\_{\text{soil}}\\ \varepsilon \sigma T^4 &= f\_{\text{c}} \, \varepsilon\_{\text{vegetation}} \sigma T\_{\text{vegetation}}^4 + (1 - f\_{\text{c}}) \, \varepsilon\_{\text{soil}} \sigma T\_{\text{soil}}^4\\ \varepsilon &= f\_{\text{c}} \, \varepsilon\_{\text{vegetation}} + (1 - f\_{\text{c}}) \, \varepsilon\_{\text{soil}} \end{aligned} \tag{25}$$

With *�* the satellite-based emissivity, *�vegetation* & *�soil* assumed fixed emissivity values for vegetation and bare soil (satellite response dependent), *σ* the Stephan-Boltzmann constant, *T* the satellite-based land surface temperature, and *Tvegetation* & *Tsoil* the pixel vegetation and soil fractions temperatures. The first equation answers to the geographical two-source proportion within the pixel, while the second equation answers to the two-source flux merging according to Yunhao et al. (2005).

Reference ET models included are Allen et al. (1998) from Cannata (2006), Priestley and Taylor (Priestley & Taylor, 1972) and Hargreaves (Hargreaves et al., 1985), Modified Hargreaves (Droogers & Allen, 2002), Hargreaves-Samani (Hargreaves & Samani, 1985). Only the reference ET from Allen et al. (1998) is being used as a precursor of SSEB (Senay et al., 2007) and METRIC (Allen et al., 2007) actual ET. It was found preponderant to have a minimum group of reference ET models available as baseline for all the work, especially when looking into geographical areas where meteorological data has always been dominant in agricultural literature.

Some models requiring operator intervention (SEBAL, METRIC) have add there internals modified with specially designed heuristics acting as operators. Initial developments were not looking into heuristics but stochastic algorithms. Some efforts using a genetic algorithm were eventually too expensive in processing time, while at the same time end-member selection information were becoming more common (Chandrapala & Wimalasuriya, 2003; Timmermans et al., 2007). Thus heuristics were designed and implemented on a regional basis, initially studied under the Greek conditions for the purpose of Alexandridis et al. (2009) and Chemin et al. (2010). Eventually, the heuristics are extended to fit data sources, continent/climate combinations and model types on an adhoc basis as new regions are included into the geographical scope of research.

### **5. Initial results**

In the case of SEBAL heuristic, the convergence reached 82% of the images processed for the Australian Murray-Darling Basin (1 Million Km2), enabling the automatic processing of 3635 MODIS multi-tiles images within a single day of computing. Fig. 3 is the output from SEBAL with such heuristic for some irrigated areas in Australia, the total area being processed amounts to more than 5 Billions pixels of ETa values, being multiplied by as many temporary rasters and original data as required for each of the ET models. The Australian irrigation system (less than 100,000 ha) has a sharp, contrasted and well-defined pattern of water depletion, characteristic of continental dry climate with high water supply control for defined periods of the year where crops are in the field.

10 Will-be-set-by-IN-TECH

skin temperature range. The available information to possibly close the equations unknown parameters are relating to satellite *T* after Price (1984) and the satellite-based emissivity (*�*), partitioned into *�vegetation* = 0.93 & *�soil* = 0.97 quoting Sui et al. (1997). Unless missing information or other sources of solution dimension constraints are published, the proposition of the model is so far considered impossible and not included in the framework until such

With *�* the satellite-based emissivity, *�vegetation* & *�soil* assumed fixed emissivity values for vegetation and bare soil (satellite response dependent), *σ* the Stephan-Boltzmann constant, *T* the satellite-based land surface temperature, and *Tvegetation* & *Tsoil* the pixel vegetation and soil fractions temperatures. The first equation answers to the geographical two-source proportion within the pixel, while the second equation answers to the two-source flux merging according

Reference ET models included are Allen et al. (1998) from Cannata (2006), Priestley and Taylor (Priestley & Taylor, 1972) and Hargreaves (Hargreaves et al., 1985), Modified Hargreaves (Droogers & Allen, 2002), Hargreaves-Samani (Hargreaves & Samani, 1985). Only the reference ET from Allen et al. (1998) is being used as a precursor of SSEB (Senay et al., 2007) and METRIC (Allen et al., 2007) actual ET. It was found preponderant to have a minimum group of reference ET models available as baseline for all the work, especially when looking into geographical areas where meteorological data has always been dominant in agricultural

Some models requiring operator intervention (SEBAL, METRIC) have add there internals modified with specially designed heuristics acting as operators. Initial developments were not looking into heuristics but stochastic algorithms. Some efforts using a genetic algorithm were eventually too expensive in processing time, while at the same time end-member selection information were becoming more common (Chandrapala & Wimalasuriya, 2003; Timmermans et al., 2007). Thus heuristics were designed and implemented on a regional basis, initially studied under the Greek conditions for the purpose of Alexandridis et al. (2009) and Chemin et al. (2010). Eventually, the heuristics are extended to fit data sources, continent/climate combinations and model types on an adhoc basis as new regions are included into the

In the case of SEBAL heuristic, the convergence reached 82% of the images processed for the Australian Murray-Darling Basin (1 Million Km2), enabling the automatic processing of 3635 MODIS multi-tiles images within a single day of computing. Fig. 3 is the output from SEBAL with such heuristic for some irrigated areas in Australia, the total area being processed amounts to more than 5 Billions pixels of ETa values, being multiplied by as many temporary rasters and original data as required for each of the ET models. The Australian irrigation system (less than 100,000 ha) has a sharp, contrasted and well-defined pattern of water depletion, characteristic of continental dry climate with high water supply control for

*vegetation* + (<sup>1</sup> <sup>−</sup> *fc*) *�soilσT*<sup>4</sup>

*soil*

(25)

*T* = *fc Tvegetation* + (1 − *fc*) *Tsoil*

*�* = *fc �vegetation* + (1 − *fc*) *�soil*

*�σT*<sup>4</sup> = *fc �vegetationσT*<sup>4</sup>

condition is being met through publication.

to Yunhao et al. (2005).

geographical scope of research.

defined periods of the year where crops are in the field.

**5. Initial results**

literature.

Fig. 3. Daily RS-based ETa (mm/day) in an Australian irrigation system (2000-2010)

Fig. 4. Daily ETa (mm/day) averages for Sri Lanka (2003-2004)

Looking into the matter of comparing ETa results from different ETa models, Fig. 4 is the averaged ETa output from two models (SEBAL and SSEB) over the tropical island of Sri Lanka in 2003 and 2004. It turns out that the relatively small island of Sri Lanka has an average ETa that is changing much more on a day to day basis than our previous example in Australia. Scale, climate, topography yield exposure to ocean events frequently, having drastic impact on thermodynamics of the island surface as the Fig. 5 also confirms. Changes between models of actual ET from SEBAL and SSEB are relatively constant throughout the RS modeling period. Actual ET from SSEB is in the upper range of SEBAL's one. The work of de Silva (1999) in the dry zone of Sri Lanka and the work of Hemakumara et al. (2003) in the wet zone of Sri Lanka are falling within the expected results found here. Likewise the average evaporative fractions found for Sri Lanka in Fig. 5 are especially leveraging the larger dry zone area of the island with value in the range of 0.3 to 0.5.

Fig. 5. Instantaneous Evaporative Fractions for Sri Lanka (2003-2004)

## **6. Conclusion**

Challenges to experimentally compare ET models are immense, the theoretical points of comparison are sometimes clear, sometimes rather difficult to pinpoint. To try and address this situation, a framework for benchmarking ET actual models has been designed. Its implementation has embedded parallel data distribution at the base of each parts of the framework to remove the resistance of the data size to process large areas, high frequency and large time period with commonly available computers.

Future work includes the finalization of SEBS (Su, 2002) and TSEB (Kustas & Norman, 1999) integration in the framework, looking for other ETa model candidates to add to existing ones. Also there is a need for designing and creating statistical tools to cross-compare several depths and layers of ETa models processing datasets. Finally, the use of OpenMPI (OpenMPI, 2011) is envisaged for concurrently running several ET models diagnostics in different multi-core machines or OpenCL (Khronos.org, 2011) kernel-based data distributed language to process all analysis as one large computation on a Graphical Processing Unit (GPU).

## **7. References**

Alexandridis, T. K., Cherif, I., Chemin, Y., Silleos, G. N., Stavrinos, E. & Zalidis, G. C. (2009). Integrated methodology for estimating water use in mediterranean agricultural areas, *Remote Sensing* 1(3): 445–465.

URL: *http://www.mdpi.com/2072-4292/1/3/445/*

Allen, R. G., Peirera, L. S., Raes, D. & Smith, M. (1998). *Crop evapotranspiration - Guidelines for computing crop water requirements - FAO Irrigation and Drainage Paper 56*, FAO - Food and Agriculture Organization of the United Nations. URL: *http://www.fao.org/docrep/x0490e/x0490e00.htm*

12 Will-be-set-by-IN-TECH

Challenges to experimentally compare ET models are immense, the theoretical points of comparison are sometimes clear, sometimes rather difficult to pinpoint. To try and address this situation, a framework for benchmarking ET actual models has been designed. Its implementation has embedded parallel data distribution at the base of each parts of the framework to remove the resistance of the data size to process large areas, high frequency

Future work includes the finalization of SEBS (Su, 2002) and TSEB (Kustas & Norman, 1999) integration in the framework, looking for other ETa model candidates to add to existing ones. Also there is a need for designing and creating statistical tools to cross-compare several depths and layers of ETa models processing datasets. Finally, the use of OpenMPI (OpenMPI, 2011) is envisaged for concurrently running several ET models diagnostics in different multi-core machines or OpenCL (Khronos.org, 2011) kernel-based data distributed language to process

Alexandridis, T. K., Cherif, I., Chemin, Y., Silleos, G. N., Stavrinos, E. & Zalidis, G. C. (2009).

Allen, R. G., Peirera, L. S., Raes, D. & Smith, M. (1998). *Crop evapotranspiration - Guidelines for*

Integrated methodology for estimating water use in mediterranean agricultural

*computing crop water requirements - FAO Irrigation and Drainage Paper 56*, FAO - Food

all analysis as one large computation on a Graphical Processing Unit (GPU).

Fig. 5. Instantaneous Evaporative Fractions for Sri Lanka (2003-2004)

and large time period with commonly available computers.

areas, *Remote Sensing* 1(3): 445–465.

URL: *http://www.mdpi.com/2072-4292/1/3/445/*

and Agriculture Organization of the United Nations. URL: *http://www.fao.org/docrep/x0490e/x0490e00.htm*

**6. Conclusion**

**7. References**

	- URL: *http://www.sciencedirect.com/science/article/B6V6C-4CYFRND-H/2/6bc33f78398c9 c74ffabe611b3ed1b6b*

URL: *http://www.springer.com/978-90-277-1247-9*

Brutsaert, W. (1999). Aspects of atmospheric boundary layer similarity under free-convective conditions, *Reviews of geophysics* 37(4): 439–451.

URL: *http://www.agu.org/journals/ABS/1999/1999RG900013.shtml*


URL: *http://www.sciencedirect.com/science/article/B6T3X-47PPC5B-1/2/d0844a2f9391840 f54c00ff0da926a75*


URL: *http://www.sab.ac.lk/journal/1999/1999A10.pdf*


URL: *http://journals.ametsoc.org/doi/abs/10.1175/1520-0450%281996%29035%3C2091% 3ARARSDS%3E2.0.CO%3B2*


URL: *http://onlinelibrary.wiley.com/doi/10.1002/hyp.7104/abstract*


14 Will-be-set-by-IN-TECH

Choudhury, B., Idso, S. & Reginato, R. (1987). Analysis of an empirical model for soil heat flux

de Silva, R. (1999). A comparison of different models of estimating actual evapotranspiration

Droogers, P. & Allen, R. G. (2002). Estimating reference evapotranspiration under inaccurate

Friedl, M. A. (1996). Relationships among remotely sensed data, surface energy balance,

Friedl, M. A. (2002). Forward and inverse modeling of land surface energy balance using surface temperature measurements, *Remote Sensing of Environment* 79(2-3): 344–354. URL: *http://www.sciencedirect.com/science/article/B6V6V-44R1BH4-K/2/57215e156bb3b7*

Gao, Y. & Long, D. (2008). Intercomparison of remote sensing-based models for estimation

García, M., Villagarcía, L., Contreras, S., Domingo, F. & Puigdefábregas, J. (2007). Comparison

Gieske, A. (2007). Numerical modeling of heat and water vapor transport through the

Hargreaves, G. H. & Samani, Z. A. (1985). Reference crop evapotranspiration from

Hargreaves, G. L., Hargreaves, G. H. & Riley, J. P. (1985). Agricultural benefits for senegal river basin, *Journal of Irrigation and Drainage Engineering* 111: 113–124. Hemakumara, H., Chandrapala, L. & Moene, A. (2003). Evapotranspiration fluxes over

Jia, L., Su, Z., van den Hurk, B., Menenti, M., Moene, A., De Bruin, H. A. R., Yrisarry, J. J. B.,

URL: *http://onlinelibrary.wiley.com/doi/10.1002/hyp.7104/abstract*

reflective and thermal data, *Sensors* 7(6): 860–883. URL: *http://www.mdpi.com/1424-8220/7/6/860/* GDAL (2011). Gdal - geospatial data abstraction library.

> Netherlands, pp. 71–83. 10.1007/978-1-4020-6218-6\_6. URL: *http://dx.doi.org/10.1007/978-1-4020-6218-6\_6*

temperature, *Applied Engineering in Agriculture* 1: 96–99.

*52f9fca7d9aef38216c*

*University Journal* 2: 87–100.

*Meteorology* 35(11): 2091–2103.

*3ARARSDS%3E2.0.CO%3B2*

*681a6460684503a761*

URL: *http://www.gdal.org*

*Water Management* 58: 109–122.

22: 4850–4869.

URL: *http://www.sab.ac.lk/journal/1999/1999A10.pdf*

URL: *http://dx.doi.org/10.1023/A:1015508322413*

data conditions, *Irrigation and Drainage Systems* 16: 33–45.

under a growing wheat crop for estimating evaporation by an infrared-temperature based energy balance equation, *Agricultural and Forest Meteorology* 39(4): 283–297. URL: *http://www.sciencedirect.com/science/article/B6V8W-4894NGY-9B/2/5668c2204c443*

from potential evapotranspiration in the dry zone of sri lanka, *Sabaragamuwa*

and area-averaged fluxes over partially vegetated land surfaces, *Journal of Applied*

URL: *http://journals.ametsoc.org/doi/abs/10.1175/1520-0450%281996%29035%3C2091%*

of evapotranspiration and accuracy assessment based on swat, *Hydrological Processes*

of three operative models for estimating the surface water deficit using aster

interfacial boundary layer into a turbulent atmosphere, *in* B. J. Geurts, H. Clercx & W. Uijttewaal (eds), *Particle-Laden Flow*, Vol. 11 of *ERCOFTAC Series*, Springer

mixed vegetation areas measured from large aperture scintillometer, *Agricultural*

Ibanez, M. & Cuesta, A. (2003). Estimation of sensible heat flux using the surface energy balance system (sebs) and atsr measurements, *Physics and Chemistry of the Earth* 28(1-3): 75–88. Applications of Quantitative Remote Sensing to Hydrology.

URL: *http://www.sciencedirect.com/science/article/B6X1W-483SMNY-5/2/735667f1cf3fef ea27e001736c7d728a*

Khronos.org (2011). Opencl: Open source computing language.

URL: *http://www.khronos.org/opencl/*


URL: *www.agu.org/journals/ABS/JD089iD05p07231.shtml*


*Hydrometeorology* 9(5): 903–919.

URL: *http://journals.ametsoc.org/doi/abs/10.1175/2008JHM920.1*


URL: *http://www.informaworld.com/10.1080/01431160512331314074*
