Spatial variability of mean daily estimates of actual evaporation from remotely sensed imagery and surface reference data

Land surface evaporation has considerable spatial variability that is not captured by point scale estimates calculated from meteorological data alone. Knowing how evaporation varies spatially remains an important issue 10 for improving parameterisations of land surface schemes and hydrological models, and various land management practices. Satellite-based and aerial remote sensing has been crucial for capturing moderate to larger scale surface variables to indirectly estimate evaporative fluxes. However, more recent advances for field research via unmanned aerial vehicles now allows for the acquisition of more highly detailed surface data. Integrating models that can estimate actual evaporation from higher resolution imagery and surface reference 15 data would be valuable to better examine potential impacts of local variations in evaporation on upscaled estimates. This study introduces a novel approach for computing a normalised ratiometric index from surface variables that can be used to obtain more realistic distributed estimates of actual evaporation. For demonstration purposes the Granger and Gray evaporation model (G-D) was applied at a rolling prairie agricultural site in central Saskatchewan, Canada. Visible and thermal images and meteorological reference data required to parameterise 20 the model were obtained at midday. Ratiometric indexes were computed for the key surface variables albedo and net radiation at midday. This allowed point observations of albedo and mean daily net radiation to be scaled across high resolution images over a large study region. Albedo and net radiation estimates were within 5 – 10 % of measured values. A daily evaporation estimate for a grassed surface was 0.5 mm (23 %) larger than eddy covariance measurements. Spatial 25 variations in key factors driving evaporation and their impacts on upscaled evaporation estimates are also discussed. The methods applied have two key advantages for estimating evaporation over previous remote sensing approaches, 1. Detailed daily estimates of actual evaporation can be directly obtained using a physically-based evaporation model, and 2. Analysis of more detailed and reliable evaporation estimates may lead to improved methods for upscaling evaporative fluxes to larger areas. 30 1 Background and introduction ‘Actual’ evaporation is the water vapour physically transferred from a surface (e.g. plants, soil or water) to the atmosphere over a given time period (e.g. hourly or daily). Reliable estimates of actual evaporation are often


20
the model were obtained at midday.
Ratiometric indexes were computed for the key surface variables albedo and net radiation at midday. This allowed point observations of albedo and mean daily net radiation to be scaled across high resolution images over a large study region. Albedo and net radiation estimates were within 5 -10 % of measured values. A daily evaporation estimate for a grassed surface was 0.5 mm (23 %) larger than eddy covariance measurements. Spatial 25 variations in key factors driving evaporation and their impacts on upscaled evaporation estimates are also discussed. The methods applied have two key advantages for estimating evaporation over previous remote sensing approaches, 1. Detailed daily estimates of actual evaporation can be directly obtained using a physically-based evaporation model, and 2. Analysis of more detailed and reliable evaporation estimates may lead to improved methods for upscaling evaporative fluxes to larger areas. 30 1

Background and introduction
'Actual' evaporation is the water vapour physically transferred from a surface (e.g. plants, soil or water) to the atmosphere over a given time period (e.g. hourly or daily). Reliable estimates of actual evaporation are often 2 needed over large spatial scales for applications such as water resource management, agriculture, ecology, and forecast modelling of weather and climate. However, estimates (or measurements) are often calculated at point scales with footprints that can range from centimetres to several kilometres or more (Brutsaert, 1982).
Consequently, point scale footprints in heterogeneous landscapes may contain large variability that needs to be considered more appropriately.

5
From an ecological standpoint, heterogeneous landscapes are comprised of distinct topographic features, land cover types, biological attributes and other physical properties that exhibit observable patterns in the order of meters (e.g. Yates et al., 2006;Zhang and Guo, 2007). Therefore, variable surface properties and state conditions exert a strong control on local surface energy fluxes. As a result, "scaling" evaporation estimates over large areas must consider a potential loss of information due to upscaling processes. The potential impacts of spatial variability 10 on larger scale estimates of evaporation is still not well understood, and previous estimation and scaling methods have not examined this issue in detail.
For example, hydrologic and atmospheric modelling applications often require large to regional scale evaporation estimates but the underlying variability is difficult to examine practically (e.g. Avvisar and Pielke (1989), Baldocchi (2005), Brutsaert (1998), Claussen (1991;, Klaassen (1992), Klaassen and Claussen 15 (1995). Courault et al. (2005) and Gowda et al. (2007) have also reviewed remote sensing approaches which integrate surface images to derive key variables needed to parameterise various energy-balance type evaporation models. This generally results in estimates at moderate scales of input images (e.g. Bisht et al., 2005). Purely empirical methods have also correlated evaporation with vegetation indices (e.g. Nagler et al., 2005).
A common remote sensing method has been to calculate evaporation indirectly as a residual of a simplified 20 energy balance (e.g. Jackson et al., 1977;Seguin et al., 1989;Bussières et al., 1996). In such cases surface temperatures derived from thermal imagery is a critical input. More complex resistance-type formulations also exist based on developments by Monteith (1965); e.g. Norman et al. (1995), Anderson et al. (1997), Boegh et al. (2002), Houborg and Soegaard (2004), and Anderson et al. (2007). However, such approaches are computationally intensive and parameterising the resistance terms is difficult without detailed data.
25 Colaizzi et al. (2006) reviewed scaling approaches based on an evaporative fraction determined through complex solar radiation modelling. Mu et al., (2007) and Fisher et al. (2008) discuss advanced methods for deriving global scale estimates based on expert knowledge and detailed data sets obtained from the Ameriflux network in association with a global Fluxnet (Baldocchi et al., 2002). As an alternative to more complex methods, Granger 3 Such analysis could then help advance understanding how the underlying surface variability may impact upscaling of point evaporation estimates to larger areas.
The goal of this study is to examine the spatial variability of key surface variables driving point scale evaporation estimates and their impact on upscaled evaporation estimates to larger areas, and includes two objectives. The first is to use the spatial variability captured from high resolution one-time-of-day visible and 5 thermal images taken at midday, to scale (i.e. distribute) energy-balance factors that drive combination evaporation models. This objective required the development of a scaling method that derives a ratiometric index of surface variables from the midday images. The resulting information can then be used to distribute a measured value of mean daily net radiation which can be used as input for deriving detailed estimates of evaporation. The second objective is to examine impacts of variations of the key surface variables on daily estimates of evaporation, 10 including point-scale and larger areal estimates.
Methods applied here integrated high spatial resolution images with an initial pixel size of less than 5 m taken during a flight on Aug 5, 2007 over a rolling prairie agricultural landscape. It also included point measurements of surface reference data, a physically-based evaporation model, and analysis using ArcGIS (v9.x) software.
Estimates of mean daily 'actual' evaporation were calculated with the complementary feedback model of Granger 15 and Gray (1989). The model is a useful alternative to more complex methods that require resistance parameters and is well suited for a variety of Canadian environments (e.g. agricultural, prairie and boreal).

Study area
A case study was conducted on Aug 5, 2007 at St. Denis National Wildlife Area (SDNWA) located in the parkland ecoregion of the Canadian Prairies in central Saskatchewan (see . Figure 1 shows a photo 20 of the study region taken during the study flight on Aug 5 to capture visible and thermal images from hand held cameras. The landscape was characterised by hummocky, gently rolling terrain, and a few slopes of up to 10 -15 %. Elevations ranged from 540 m to 565 m and land use consisted of mixed cool season grasses, tall brome grass, cultivated land, and wetlands surrounded by trees or dense grass; all grasses and crops were C3 types. The soil region is classified as dark brown chernozem and soil texture is predominately silty loam (van der Kamp et al., . Granger and Gray (1989) developed the 'G-D' model from the complementary relationship of Bouchet (1963) and Penman's (1948) combination model. The G-D model extends the potential evaporation model to non-saturated 30 surfaces using the relative evaporation term, G. This is defined as a ratio of actual to potential evaporation which depends on the relative drying power of the air, D. The underlying theory is based on the reduction of water availability as a surface dries but the 'potential' evaporation increases due to a subsequent rise in surface temperature.

4
The development of G eliminated the need for observations of surface temperature and vapour pressure. As a result, estimates of actual evaporation were obtainable for non-saturated surfaces with atmospheric data alone (Granger, 1989). The general form of the equation is: (1) The available energy term is driven by net radiation, Q* (W m -2 ) calculated as a sum of the net shortwave and 5 longwave radiation components and the ground heat flux, Qg (W m -2 ), and slope of the saturation vapour pressure curve, Δ.
The aerodynamic term includes the psychrometric constant, γ, and "drying power of the air", EA, calculated using a Dalton type formula: 10 where f(u) is a vapour transfer function, and the atmospheric vapour pressure deficit at 2 m height is derived from e*a (saturated) and ea (actual). Pomeroy et al. (1997) empirically derived f(u) as a function of wind speed and aerodynamic roughness length, z0 (m) from extensive field data collected for prairie, boreal forest, and northern cold region environments in western Canada:

15
where u is the mean daily wind speed (m s -1 ).
The G-D complementary feedback method is driven by the non-linear relationship between G and the relative drying power of the air, D: where D is a function of the humidity deficit and available energy given by: and λ is the latent heat of vaporisation (kJ kg -1 ).
The feedback model approach is generally applicable for regions indicated above and can be used when detailed soil moisture information is lacking, except under conditions of severe moisture stress.  evaluated point scale evaporation estimates obtained with the G-D model over mixed grasses at the same study 25 area under conditions of non-limiting soil moisture using field data collected during 2006. Estimates obtained with the G-D model compared well with eddy covariance (EC) measurements and was less data intensive than the Penman-Monteith model which also performed well in that study. Armstrong et al. (2010), however, found the use of appropriate soil moisture constraints was required to obtain reliable estimates with the G-D model when applied under drought conditions to a mixed prairie site at Lethbridge, 30 Alberta, Canada.

Field observations
Inputs needed to parameterise Eq. (1-5) included mean daily net radiation, air temperature, humidity, wind speed, and surface roughness length. Outgoing radiation components were derived from high resolution digital (Canon Powershot A70 camera) and thermal (FLIR ThermaCAM P20 imaging radiometer) images taken on Aug 5, 2007 during a Cessna flight at midday. Meteorological and surface observations (radiation components, temperature, 5 humidity, wind speed) were recorded at two station locations. Figure 1 shows the relative locations which included one reference site and one validation site. An independent station operated by the National Water Research Institute (Environment Canada, Saskatoon) was also available which provided a second validation site for estimates of daily net radiation. Air temperature, humidity, wind speed, incoming radiation and surface temperature were recorded  Disney et al. (2004) and Zhang and Guo (2007). Samples were taken at 4.5 m intervals along a site transect (see

Surface reference meteorological parameters
At midday, reference observations of the shortwave, K↓ (835 W m -2 ) and longwave, L↓ (320 W m -2 ) irradiance and albedo (0.153 over grass) were obtained from the CNR1. Both K↓ and L↓ were assumed to be uniform over the field given the images used for calculations were cloud free. The reference values were used in the calculation of midday normalised ratiometric indexes for albedo and net radiation. A mean daily reference value of 155 W m -35 2 was obtained for Q* which was the input for distributing the final estimates of the mean daily net radiation using the ratiometric index (discussed in Section 3.4).

6
The daily ground heat flux measured at two locations was relatively small, 8 W m -2 and 2 W m -2 , due to the continuous and tall grass canopy cover. For this study the ground heat flux was ignored, in part due to the small measured values, but also to limit errors introduced from using a numerical solution to estimate ground heat flux, which would require more detailed cover information at each image pixel. To standardise parameterisation of the G-D model for calculating mean daily estimates of evaporation, the reference values of mean daily air temperature 5 (19.6 °C), humidity deficit (1.1 kPa) and wind speed (3 m s -1 ) were also taken to be uniform. Differences in mean daily air temperature between the measurement locations were approximately 1 °C.
Potential impacts of assuming uniform humidity and wind speed for parameterising the aerodynamic term was also considered. This was done by examining G-D evaporation estimates derived from 2006 field data collected over different surfaces at the same study area . Implications of neglecting the ground heat flux and treatment of the meteorological parameters (above) with respect to modelling uncertainty are given in the discussion.

3.4
Deriving key surface variables from one-time-of-day images 20

Theoretical basis for a normalised ratiometric index for surface properties
Remotely sensed images contain valuable information that characterise highly variable surface properties (e.g. reflectance, temperature, RGB and greyscale DNs, etc.). Theoretically, relative variations in these properties can be quantified using a "normalised" ratiometric index. For example, consider the 'evaporation ratio' (ER) defined simply as the ratios of individual evaporation rates at different spatial locations (Ei) to a reference rate obtained at 25 a specified location (Eref): At the reference location ER = 1 and will be the same at any other locations where Ei = Eref, but will vary from unity at all other spatial locations.
For obvious reasons the theory integrates well for computations with pixel-based images. Subsequently, an 30 evaporation rate measured at the reference pixel could be scaled to all other pixels by multiplication with the value of ER at each pixel. The concept may be extended to other surface variables (e.g. albedo, net radiation etc.) required to parameterise an evaporation model. 7 3.4.2 Distributing surface variables using a normalised index of relative surface ratios Methods described here assume spatial variations in surface variables driving net radiation (and driving evaporation) are near their maximum around solar noon. This is likely to be valid within 2 hours from the actual time of solar noon (Colaizzi et al., 2006). Net radiation is the major component needed to determine available energy for the conversion of water to an equivalent depth of water vapour. Net radiation is determined from the 5 shortwave and longwave radiation components measured in the electromagnetic spectrum between approximately 0.3 -4 μm and 4 to 14 μm respectively (Liang, 2004;Zoran and Stefan, 2006).
Traditionally, radiative terms needed for estimating evaporation are derived from satellite-based imagery, e.g.
Landsat, AVHRR (Advanced Very High Resolution Radiometer) and MODIS (Moderate-resolution Imaging Spectroradiometer). However, satellite-based methods are often limited by cloud contamination, varying spatial 10 and temporal resolutions and sensor footprint mismatches. Under clear skies it can be assumed shortwave and longwave irradiance are uniform over an agricultural field area. Very high resolution images taken near the surface could be used to derive the much more variable surface reflected shortwave and emitted longwave radiation components.

15
can be calculated at every pixel location within visible and thermal images using: where subscript 'i' is an individual value at each pixel and 'ref' is the value at the reference pixel. Incoming 20 shortwave (K↓) and longwave radiation (L↓) components can reasonably be assumed to be uniform over the field under clear skies, so αR and L↑R can be further integrated to derive a ratiometric index of the midday net radiation, Subsequently, single measured values of albedo and mean daily net radiation taken at a reference pixel can be 25 scaled (i.e. distributed) across all other pixels via multiplication with the surface ratios derived for Q*R. The next section illustrates the indexing method for deriving accurate estimates of albedo.

Ratiometric index method for albedo estimates from digital visible images
Surface albedo (α) represents a crucial radiation loss term for radiative transfer calculations and surfaceatmosphere energy and mass exchanges (Sellers et al., 1997;Liang, 2000;Lucht et al., 2000;Roberts, 2001;Liang 30 et al., 2003;Disney et al., 2004;Liang, 2004. Its calculation can be complex and is typically a major source of uncertainty (Yang et al., 2008).

8
In this case, a measured value of broadband albedo obtained at a reference pixel location could be scaled to every other pixel within a high resolution visible image. In 2007, digital photos were taken with a Canon Powershot A70 camera; max resolution 2048 x 1536 pixels, CCD (charge-coupled device) imager, DIGIC (Digital Imaging Core) processor. This resulted in very high resolution (< 1 m pixels) 'visible images' without cloud cover issues.
While digital cameras may not cover the full visible and near-infrared spectrum of advanced measurement 5 sensors, the imaging techniques are based on the same principles. Corripio (2004) demonstrated how accurate estimates of snow albedo could be obtained from digital images using a linear scaling technique applied to a measured albedo at a reference pixel location. A key step for our analysis was to transform an RGB digital photo to a single band 8-bit greyscale image with DNs (Digital Numbers) ranging from 0 -255.
This resulted in higher DNs associated with more reflective surfaces (i.e. brighter) and lower DNs with less 10 reflective surfaces (i.e. darker) in accordance with principles for calculating albedo. The resulting albedo map was aggregated to a pixel size of 5 m for a practical analysis and georectified to 100 GPS ground control points. The study area included some dried wetlands and others containing open water. This allowed for application of a simple Dark Object Subtraction (DOS) method to correct potential atmospheric effects (Song et al., 2001;Liang, 2004).
Applying Eq. (11) made it possible to accurately estimate albedo at individual pixel locations, αi, from a 15 measured broadband value taken at the reference location, αref, and a ratiometric index for DNs calculated at each pixel: where DNi is the digital number of an individual pixel and DNref is the reference pixel value.
Angular measurements of broadband albedo from 0.3 -3.0 μm with a hemispherical CNR1 directly accounted 20 for bidirectional reflectance properties of the mixed grassed surface. Therefore, the point observations satisfied the bidirectional reflectance considerations related to albedo estimation techniques (Nicodemus et al., 1977;Lucht et al., 2000;Roberts, 2001). Figure 2 shows a sample of measured reflectance values from a grassed upland area. The spectral reflectance from the mixed grasses was virtually identical to the response for healthy (green) winter wheat (see Disney et al., 25 2004). For wheat, they concluded reflectance was directionally invariant and reasonably could be assumed to behave as a Lambertian surface (i.e. scattering light equally in all directions). The field measured spectral reflectance over the grassed surface was therefore treated similar to a Landsat measured reflectance and divided into respective wavelengths for Landsat wavebands 1, 3, 4, 5, and 7.
An empirical linear approximation for narrow-to-broadband albedo conversion applicable to Landsat imagery 30 (see Liang, 2000) was then applied to the field measured spectral reflectance data. This allowed for a direct comparison of albedo estimates and measurements along the sampling transect shown in Fig. 3. Due to the 4.5 m sampling distance some pixels contained a single value whilst others contained two values, in which case the mean value was used for comparison. 9 3.4.4 Surface temperature (emitted longwave radiation) For the case of surface temperatures the indexing method was not needed because high resolution observations were obtained directly with a hand held Forward Looking Infrared (FLIR) ThermaCAM P20 imaging radiometer.
The P20 used a Focal Plane Array, uncooled microbolometer, with a maximum image resolution of 320 x 240 pixels, a 24° by 18° field of view, and spatial resolution of 1.3 milliradians. The spectral range was 7.5 -13 μm 5 which is similar to traditional satellite sensors (e.g. Landsat, MODIS, and AVHRR 10 -12.5 μm, and ASTER 8 -12 μm).
A standard emissivity of 0.98 was assumed for the cover types encountered, and ambient air temperature / humidity and distance between the surface and camera detector were set based on observations. The Stefan-Boltzmann equation was applied to transform surface temperatures into values of outgoing longwave radiation 10 needed to estimate the net radiation. For the flight height of approximately 1 km the FLIR produced a surface pixel resolution of approximately 3 m. The longwave radiation map was then aggregated to 5 m resolution and georectified using the resulting map of albedo estimates as the reference.

Surface roughness length
The study area is characterised by a complex landscape with rolling terrain covered by tall grasses, cultivated land 15 and narrow rings of tall grasses, tall shrubs or tall trees surrounding ponds and dried wetlands (see Fig. 1 the standard deviation and associated DN for each pixel is then sorted (low to high) and a bin range assigned; a 35 class width tolerance was set for pixels having similar standard deviations and all values within a specified range were assigned to the same class. Where pixel values were outside the range, but class boundaries overlapped, a mid-point was determined and a new class created.
The initial 13 classes were manually reclassified into three general classes (fallowed/cropped, grass, and trees) based on the extents of the dominant land cover types observed in the original and classified images. Characteristic roughness lengths, zo were then selected for each class based on values reported for similar surface types (Brutsaert, 5 1982). In this case, a value of 0.05 m was used for the fallowed/cropped class and 0.10 m for the taller more dense grass class. Narrow rings of shrubs and trees around wetlands were also dense and much taller than the surrounding cover with heights varying from 3 m to 10 m. Brutsaert (1982) reported roughness lengths between 0.2 and 0.4 m for similar type surface elements. The value of 0.40 m was chosen to reflect an expected increase in roughness and turbulent exchange due to a likely reduction in wind speed compared to the surrounding cover types.

Results and discussion
Calculating detailed daily evaporation estimates and examining spatial scaling issues required the midday inputs and temporal scaling function to be computed first. Analysis of midday inputs are described in sections 4.1 -4.4.
Sections 4.5 and 4.6 discuss the sensitivity of the G-D model to the midday evaporation ratio and the temporal 20 transfer function required for scaling a single measured value of mean daily radiation across a midday image.
Sections 4.7 -4.9 discuss the accuracy of the resulting evaporation estimates, variations in statistical distributions of key driving components, and scaling implications for larger scale evaporation estimates.

Validation of albedo estimates
At midday on Aug 5, 2007 a reference albedo, αref = 0.153 was obtained from observed irradiance and reflectance 25 of shortwave radiation over the mix of green grasses. The albedo map resulting from Eq. (11) and locations of reference and validation points is shown in Fig. 3. Vegetation was similar at both EC station locations and the scaled albedo estimate (0.164) agreed well with the measured value (0.167) at the validation site.
Validation of mean and range of albedo estimates obtained for major land covers in Fig. 3

is summarised in
Tab. 1. Estimates of albedo from the image compared well with values expected for grasses, agricultural crops, 30 deciduous trees, and bare soils reported in Brutsaert (1982). The root mean square error (RMSE) of albedo estimates from the pixels and measured albedo values was approximately 3.5 % which is within an expected error of 2 % to 5 % for research purposes (Liang, 2004).
The 5 m resolution albedo map highlighted surface information within the landscape that would be more generalised using coarser data. For example, the albedo map depicted distinct boundaries separating regions dominated by brome (BG) and mixed grasses (MG), cultivated/crop area (C), sparsely vegetated, fallow area (F), and wetland fringe vegetation (WL). Also, wetland extents and fringe vegetation were observable in areas surrounded by other vegetation types (see Fig. 1 for reference). The detailed variations are expected to impact 5 evaporation estimates through relative increases at pixels with relatively lower values of albedo which results in an increase in the available energy, whereas relatively higher values will reduce the available energy and estimates of evaporation. Figure 4 shows the resulting map of emitted longwave radiation, L↑ derived from the image of observed surface 10 temperatures, with estimates ranging between 380 W m -2 to 480 W m -2 . At the two EC station locations a comparison was made among the midday surface longwave measurements from the P20 camera, a narrow-beam Exergen infrared thermocouple (IRTC) radiometer and surface emitted longwave obtained from the CNR1.

Validation of longwave radiation estimates
The P20 measurement compared well with the IRTC values with differences less than -12 W m -2 . Compared with the CNR1 values the differences where slightly larger at -30 W m -2 but was still within 8 % error. Generally, 15 the differences were small considering relative changes in midday surface net radiation and the magnitude of incoming radiation components is considerably larger. The variable footprints of the different measurements and absorption properties of water vapour which can reduce the signal are likely sources of error. Also, dust and heating of the CNR1 downward facing pyrgeometer may introduce another source of error.
For estimating evaporation, the smaller values of emitted longwave resulting from lower surface temperatures 20 will increase the available energy relative to larger values which imply reduced water availability and evaporative cooling.

30
Normalised indexes for albedo, surface temperature and roughness were considered to examine the relationships of variations in these variables compared to the evaporation ratio. Figure 6 shows the expected physical behaviour of these variables within the G-D model. Only the actual range of values computed for the ratiometric index was considered so the physical variations could be shown more clearly. In this case, the impacts of relative changes in the apparent inverse linear relationships between ER and αR and L↑R, and slight non-linear relationship between ER and zoR on evaporation estimates can be computed.
For example, a relative increase in the surface temperature via L↑R of 0.18 (or 18 %) was shown to reduce ER by 10 %. By comparison, an increase in reflected shortwave radiation via αR of 0.30 (or 30 %) reduced ER by only 5 %. In the case of surface roughness, an increase of 250 % was needed to reduce ER by 10 %. Consequently, a 5 relative reduction in evaporation rates may be expected where albedo, surface temperature and surface roughness tends to be larger. In the latter case, step changes in surface roughness would increase the relative drying power, D but the relative evaporation, G results in a relative decrease due to the inverse non-linear relationship. The increased sensitivity of ER to L↑R also indicates detailed spatial variations in surface longwave radiation is an important factor for estimating evaporation.

Temporal transfer function: normalised ratiometric radiation index
Net radiation is known to vary dynamically on a sub-daily basis. Eq. (10) shows how a radiation ratio, Q*R can be used as a temporal transfer function to scale estimates of mean daily net radiation over an image. In order to scale the normalised net radiation ratios from a temporal "point" at midday to a mean daily value it was necessary to examine whether a stable proportionality existed between measured values at midday and mean daily net radiation.

15
Verification of a proportionality under clear skies would eliminate the need for an empirical scaling function.
Historical records were examined for a period from May 1 -Sept 1 at three Canadian prairie locations at similar latitudes (49° -52.2°). The analysis included two field seasons at the SDNWA study site (2006)(2007). Archived data was also obtained at two short grass prairie locations; an Ameriflux network site at Lethbridge, AB (1999AB ( -2004, and Kernen Farm located at Saskatoon, SK (1999SK ( -2000. Figure 7 shows a moderately strong relationship 20 at each location; r 2 = 0.54 -0.6. These results suggested the relationship between midday and daily net radiation was stronger when midday net radiation exceeds 400 W m -2 , which is more likely for generally cloud-free days.
Confirmation of an existing proportionality suggests daily net radiation, Q*d could be scaled to all pixels across an image from a single measured reference value, Q*dref, as a function of the normalised ratio index for Q*R at midday, stated here as: Q*dref was assigned a daily measured reference value of 155 W m -2 obtained from the CNR1 at the reference station. Q*R was derived from Eq. (10) using the measured reference values of incoming shortwave and longwave radiation components, as well as the albedo and emitted longwave maps as input. The resulting mean daily net radiation map and a comparison of estimated and measured Q*d at two validation locations is shown in Fig. 8

Calculating direct estimates of actual evaporation with the G-D model
The mean daily net radiation map derived from Eq. (12) was supplied as input to the G-D model to estimate the mean daily "actual" evaporation at every land surface pixel. The reference humidity and wind speed values 5 (discussed earlier), and the surface roughness length map were used for calculating the aerodynamic terms in Eq.
(2 and 3). An estimate of mean daily evaporation was then calculated at every pixel location based on Eq. (1).
The resulting map of distributed actual evaporation estimates is shown in Fig. 9. Due to the level of detail captured in the net radiation map of Fig. 8, the map shows a physically realistic pattern of evaporation for the range of land cover types. For example, areas where vegetation was less dense (e.g. fallowed/cropped) and soil 10 surface conditions were drier showed lower rates of mean daily evaporation and higher rates were associated with the more dense grassed areas. The highest evaporation estimates were obtained where water availability is expected to be higher; e.g. among wetland fringes and depressions where albedo and surface temperatures were lower, which can be attributed to increased water availability and evaporative cooling.
The prevailing wind direction for the day (Udir) was from a north-northwest direction. A mean evaporation 15 estimate of 2.7 mm was obtained along a linear transect for pixels containing a brome grass surface immediately upwind of the EC station (Fig. 9.). The estimate of 2.7 mm was 0.5 mm or 23 % higher than the EC measured mean evaporative flux of 2.2 mm/day. Based on the cumulative flux model of Schuepp et al. (1990) it is expected that 80% of the cumulative flux would come from within an upwind distance of 100 m from the EC station along a similar path to the linear transect. The estimate of 2.7 mm/day might be lower if ground heat flux was found to 20 be a factor for reducing the available energy. However this would need to be reliably accounted for at every 5 m pixel.
On the study day the mean daily energy fluxes in W m -2 (including latent, LE and sensible heat, HE) and ratio of the energy balance closure at the validation site was: (LE + HE) / (Q* -Qg) = (63 + 55) / (144 -2). This results in an energy balance closure of approximately 83%. The resulting Bowen Ratio of 0.87 was reasonable for 25 the drying conditions and later timing in the growing season for the grasses in this semi-arid region. In this case it is expected there will be uncertainty in both estimated and measured EC fluxes due to modelling and measurement errors and different footprint scales. For instance, it is possible that LE and HE could be under-measured or the ground heat flux more variable upwind of the EC station than estimated.
Measured evaporation rates were not available for trees (dominated by Aspen) in 2007 but the G-D estimates 30 were compared against archived values from BOREAS data for an Old Aspen site from August, 1996. The G-D mean daily values in the order of 3 mm were reasonable compared to evaporation from similar trees reported by Hogg et al. (2000).
In general, the results are instructive because the G-D method has been demonstrated to provide a reasonable estimate of evaporation from the remote sensing images taken over a complex landscape.  35 have also shown similar accuracy for daily and multi-day periods from meteorological data alone during the 2006 field season at the same study area. More importantly, the ratiometric index method of scaling may be valid for 14 distributing surface radiation components across remote sensing imagery obtained from satellite, planes or near surface aerial platforms such as UAVs.

Distributions of evaporation and driving surface variables
The following sections briefly discuss the statistical distributions of the driving variables obtained from the images used for estimating evaporation and their impacts on the resulting evaporation estimates. The benefit here is the 5 statistical distributions of evaporation and key driving factors are considered to be physically meaningful. Figure 10 shows the frequency distributions of evaporation estimates and relative contributions for the energy balance and aerodynamic components. G-D model evaporation estimates appear normally distributed with a mean of 2.8 mm/day and a relatively low coefficient of variation (cv = 0.07). Distributions for the energy and aerodynamic terms were notably different due to interactions among the driving variables and different roughness 10 classes. For the energy component the distribution was continuous with a mean of 1.1 mm and cv of 0.11. The distribution mean for the aerodynamic term was 1.7 mm and was bimodal due to the differences in magnitude of the discrete roughness classes used. Figure 11 shows boxplot summaries of distributions for albedo (α), surface temperature (Ts), net radiation (Q*), relative evaporation (G), mean daily evaporation (E), and G plotted against net radiation. In the case of α,

15
values below approximately 0.10 and lower "outliers" can be attributed to uncertainty associated with the wetland vegetation and likely where surface water was not completely masked out near the edges of ponds. In general, the distributions for α and Ts which drive the net radiation appear to show more notable variability compared to Q* and the resulting estimates of E. For example, the distributions of α and Ts showed opposite skewness, -0.83 and 0.36 respectively. The variability for α (cv = 0.19) was larger than Ts (cv = 0.14), which is also much larger than 20 apparent variability of both Q* (cv = 0.65) and E (cv = 0.66). Mean values of the distributions for α, Ts, Q* and E were close to the median values shown within the boxes depicting the interquartile ranges.
The resulting boxplot for G and plot of G against net radiation reflects the interaction across the discrete roughness length classes. The statistical distribution of G was not continuous due to the larger step change in roughness length from 10 cm to 40 cm. Plots of G against net radiation resulted in three distinct linear relationships 25 which can be attributed to using a uniform mean daily wind speed and humidity deficit for calculating the drying power of the air, EA. In this case, the potential variability of G associated with the relative drying power, D is limited to variations in roughness length, but may vary more with changes in wind speed and the humidity deficit.

Spatial variations within roughness length classes
Boxplots shown in Fig. 12 further characterise the spatial variability of α, Ts, Q*, G and E within the three roughness length classes. In general, the plots depict notable shifts in the interquartile ranges of the distributions across the classes. Results of paired Kolmogorov-Smirnov tests for each variable showed statistical differences among the distributions for each roughness class where highly significant (p-value < 0.001). General increases in 35 albedo and surface temperature are clearly shown across the 40 cm to 20 cm and 10 cm roughness length classes.
Consequently, there was a notable reduction in net radiation and similar reduction in evaporation estimates moving from the higher to lower roughness length classes. Figure 12 also clearly shows evidence of a non-linear, inverse relationship between G and net radiation across the roughness length classes. In this case, the interaction between G and Q* appears to offset the potential increase in mean estimates of E associated with increases in available energy. In the other words, the impact here is more likely to be a reduction in the variability of the 5 evaporation estimates.

Scaling implications
There is a potential for areal estimates of evaporation to vary depending on how upscaled estimates are calculated For the entire area and also each roughness length class, mean values of the driving factors were derived and evaporation estimates were recalculated using Eq. (14). The relative evaporative contributions attributed to the energy and aerodynamic terms are provided in Tab. 2. For the different roughness length classes the range of 35 16 evaporation estimates attributed to E_energy was only 0.2 mm/day and nearly 0.7 mm/day for E_aero, and the difference in total evaporation, E_total was 0.5 mm/day. A bias toward larger evaporation estimates might be expected given the increase in energy availability and enhanced turbulence with an increase in roughness.
However, the potential bias was offset by the interaction of G and Q*, and also G and EA. Table 2 also compares evaporation estimates calculated from only mean values of G and Q* and EA to the 5 expected rates for each roughness length class. "Expected" rates were calculated from the mean values for all pixels assigned within each roughness length class. Evaporation rates derived from the mean input values alone were only between 0.14 mm and 0.2 mm less than the expected mean rates for each roughness length class.
Upscaling the driving factors to the entire area also had no impact on the resulting estimate as the difference was just 0.1 mm. 10 4.9 Examining spatial covariance among key variables Whether evaporation estimates might be influenced by a spatial covariance between driving factors was also examined. The Pearson correlation coefficient, r was used to evaluate correlations among the driving factors distributed over the field area. By definition, Pearson's correlation is the ratio of the covariance (the numerator) between two variables normalised by the product of their standard deviations as follows: where Xi and Yi are the respective values of the variables, the overbar denotes the mean value and n is the number of pairs. A strong correlation between two variables might suggest the existence of covariance that could influence upscaled estimates of evaporation. Given the roughness classes used represent discrete data, and a lack of more detailed meteorological data to parameterise the aerodynamic term, further evaluation related to climate factors 20 would not be meaningful. So only an examination of the factors driving the energy term are considered here.
In this case, a potential covariance between G (dimensionless) and Q* (expressed in mm/day) can be considered directly. By rearranging Eq. (15) the covariance can be obtained by multiplying the correlation coefficient and the product of the standard deviations of Q* (0.34 mm/day) and G (0.021 mm/day). The correlation between Q* and G over the field area produced a coefficient, r = -0.67.

25
When multiplied in series (r = -0.67)*0.34*0.02 this results in a covariance of approximately -0.0049 mm/day. This suggests interactions between Q* and G for the methods applied here would have no further statistical influence on upscaled evaporation estimates. Unfortunately, no comment can be made regarding covariance between G and the turbulent flux component in the G-D model. Such analysis would require more detailed observations of air temperature, humidity and wind speed, and possibly a more sensitive combination model.

30
Nevertheless, the impact of combined interactions within the G-D model was to effectively reduce the overall variability of point scale evaporation estimates and also upscaled estimates derived from different computation methods for parameterising the model.

General uncertainty of methods applied
Implications of general modelling assumptions and uncertainty of the methods applied are briefly discussed here.
The daily ground heat flux was considered to be negligible for the study day as the mean daily flux was relatively small at the two stations where there was good canopy coverage. Where local cover properties are similar (e.g. brome grass and mixed grass areas) the ground heat flux may be of a similar magnitude but would be greater where 5 surface cover is less dense or where water availability is limited, due to differences in the surface properties.
Further, the relative evaporation, G which is based on the relative drying power, D integrates available energy which in this case directly includes the spatial variability of surface temperature. As such, evaporation might be overestimated in areas with more ground exposer and higher surface temperatures, where ground heat flux may be appreciably larger (e.g. fallowed/crop area). This interaction could result in a further reduction in the energy 10 available for the evaporation estimate.

Summary and conclusions
This study examined spatial associations and physical interactions amongst key surface variables driving actual evaporation estimates, and impacts of their variations on various methods of upscaling estimates to a larger area.

20
The methods applied demonstrate how measured reference values of albedo and mean daily net radiation can be scaled accurately across a large field area for the purpose of deriving point scale estimates of evaporation at each pixel. This was achieved by computing a normalised ratiometric index of surface radiation from highly detailed midday visible and thermal images. At two validation sites estimates of daily net radiation showed good agreement with measured values to within 4 % and 8 % error.

25
Estimates of mean daily actual evaporation were calculated at 5 m resolution with the G-D model. The daily evaporation rate estimated for a transect upwind of the EC station was 2.7 mm for Aug 5, 2007 which was 23 % larger than the EC measured flux of 2.2 mm. Offsetting interactions between the relative evaporation term and key surface variables effectively reduce the spatial variability of evaporation estimates. As a result, differences amongst computed areal evaporation estimates were relatively small regardless of the different methods used to 30 derive representative average values of driving factors to parameterise the model, or an average evaporation rate derived from detailed estimates across the field. There was no evidence of a spatial covariance between the spatial distributions of net radiation and G, so no correction factor was identified for improving upscaled evaporation estimates.
The scaling methods applied to the energy terms using ratiometric indexes derived from detailed images could 35 generate useful diagnostic information at other study locations and potentially over much larger areas. The methods applied here may also be instructive toward improving techniques for upscaling evaporation estimates to larger areas via traditional remote sensing or climate modelling, or using a different combination model. It is expected the methods can be applied to visible and thermal images taken from cameras and sensors on a variety of sensing platforms (e.g. Satellites, planes and UAVs).       E_total is the combined total. The mean value of the distributed estimates is given by "Expected" and the difference between the total and expected is given by "Diff".