Hydrology and Earth System Sciences Global Land-surface Evaporation Estimated from Satellite-based Observations

This paper outlines a new strategy to derive evaporation from satellite observations. The approach uses a variety of satellite-sensor products to estimate daily evaporation at a global scale and 0.25 degree spatial resolution. Central to this methodology is the use of the Priestley and Taylor (PT) evaporation model. The minimalistic PT equation combines a small number of inputs, the majority of which can be detected from space. This reduces the number of variables that need to be modelled. Key distinguishing features of the approach are the use of microwave-derived soil moisture, land surface temperature and vegetation density, as well as the detailed estimation of rainfall interception loss. The modelled evaporation is validated against one year of eddy covariance measurements from 43 stations. The estimated annual totals correlate well with the stations' annual cumulative evaporation (R = 0.80, N = 43) and present a low average bias (−5%). The validation of the daily time series at each individual station shows good model performance in all vegetation types and climate conditions with an average correlation coefficient of R = 0.83, still lower than the R = 0.90 found in the validation of the monthly time series. The first global map of annual evaporation developed through this methodology is also presented.


Introduction
Detecting changes in the hydrological cycle is essential if we are to predict the impacts of climate change.However, climate change is acting on a dynamic three-dimensional globe where changes in one region may produce impacts in another.Therefore, there is a need to expand the current climate change studies to encompass the entire globe.
Correspondence to: D. G. Miralles (diego.miralles@falw.vu.nl)Precipitation and evaporation are the two key components of the global water cycle.Evaporation can cause feedbacks on large scale water processes (e.g.Poveda and Mesa, 1997) and affect the dynamics of the atmosphere due to changes in the Bowen ratio (e.g.Dow and DeWalle, 2000).While our capability of observing precipitation has considerably improved with the deployment of dedicated satellites such as the Tropical Rainfall Measuring Mission (TRMM) and in the near future the Global Precipitation Measurement (GPM), our capability of observing the return-flow of moisture from the land to the atmosphere is still poor (Dolman and De Jeu, 2010).Model estimates put the amount of evaporation from the global land masses somewhere between 58-85 10 3 km 3 yr −1 , although the exact magnitude and spatial and temporal variability are still highly uncertain (Dirmeyer et al., 2006).
If we are to effectively manage adaptation to climate change, the uncertainty in predictions of future climate must be reduced.This creates the need for evaporation products that can be used to validate components of Global Circulation Models (GCM) and serve as an observational benchmark for GCM developers (Blyth et al., 2009).The development of evaporation data sets from hydrological models, land surface parameterisation schemes, and/or through the application of the currently available data products (including remote sensing data) are therefore essential to improve predictions of future climate.
In the last two decades several attempts have been made to build global evaporation products based on a range of approaches tailored to specific input data.They can be categorized in four groups depending on whether they are based on: (1) off line models (e.g.GSWP -Dirmeyer et al., 2006), (2) remote sensing observations (Fisher et al., 2008;Jiménez et al., 2009), (3) reanalyses (e.g.ERA-Interim - Simmons et al., 2006), or (4) upscaling of in situ observations (MTE -Jung et al., 2009).Few of the existing approaches have been adapted to the global scale and daily frequency and Published by Copernicus Publications on behalf of the European Geosciences Union.D. G. Miralles et al.: Global land-surface evaporation have their results publicly available.The majority of them lack any emphasis on estimating rainfall interception loss or do not couple transpiration with observed soil moisture conditions; only a few of them (e.g.Fisher et al., 2008) include observation-based moisture constraints within their scheme.Acknowledging the differences between these approaches, in 2008 the LandFlux initiative of the GEWEX Radiation Panel raised the importance of evaluation and inter-comparison of the existing land evaporation products (Jiménez et al., 2009) towards the creation of reliable evaporation benchmarks.
The present paper outlines a new methodology to estimate global land-surface evaporation mainly based on satellite observations.The approach relies on the potential of the existing satellite-based data sets conferred by their observational nature (as opposed to modelled fields) and their ability to provide global spatial estimates (as opposed to in situ observations).The ultimate goal is to derive a global, 24 year, 0.25 degree, daily data set that can be used for studies of the global hydrological cycle.Central to the approach is the use of the Priestley and Taylor (PT) (1972) evaporation model.Because the PT equation requires a small number of inputs, and the majority can be directly observed by satellites, this strategy minimizes the number of modelled variables.Key distinguishing features are the use of microwave-derived soil moisture, land surface temperature and vegetation density, and the detailed estimation of rainfall interception loss.

Methodology
The model, known as GLEAM (Global Land surface Evaporation: the Amsterdam Methodology), is designed to maximize the use of satellite-derived observations to create a spatially coherent estimate of the evaporative flux over land.For this reason, parameterisations are chosen that have global validity; whenever possible, constant parameters are preferred over those which vary across the globe.As a consequence, the methodology distinguishes only three sources of evaporation based on the land surface type: (1) bare soil, (2) short vegetation, and (3) vegetation with a tall canopy.The snow and ice sublimation is estimated for the pixels covered in snow through a separate routine.The contribution of lakes and rivers is not modelled; the predicted evaporation therefore refers only to the land fraction of the total surface area of each grid cell.The land evaporation (E) of each gridbox is the sum of the evaporation modelled for each of the three land surface types (s), weighted by their fractional coverage (a): (1) The global model is composed of four modules.In the first module, the evaporation of intercepted rainfall from forest canopies is calculated.A separate module describes the water budget that distributes the incoming precipitation (rain and snow) over the root-zone.In a third module, the stress conditions are parameterised as a function of the root-zone available water and dynamic vegetation information.Finally, the evaporation from each of the three surface components is calculated based on the PT equation, the modelled stress, rainfall interception and snow sublimation.Figure 1 gives an overview of the structure of GLEAM and its main inputs and outputs.The interception model has already been described and validated by Miralles et al. (2010).The entire evaporation methodology is validated in the present paper.

Canopy interception loss
The neglect of evaporation from wet forest canopies, referred to as rainfall interception loss, is thought to be one of the main factors contributing to the uncertainty of global evaporation estimates (see Jiménez et al., 2011).In GLEAM, it is explicitly modelled according to Gash's analytical model (Gash, 1979;Valente et al., 1997).Following this approach, the volume of water that evaporates from the canopy is derived from the daily rainfall using parameters that describe the canopy cover, canopy storage, and mean rainfall and evaporation rate during saturated canopy conditions.A novelty in this application is the use of a remotely-sensed lightning frequency product to define global maps of monthly climatology of rainfall rate.The derivation of the parameters, validation and global implementation of the GLEAM interception model is fully described by Miralles et al. (2010).The study showed a strong correlation (R = 0.86) and a negligible bias between modelled and observed values of interception as reported in 42 field studies over different forest ecosystems.

Soil water content
The second module computes a daily running water balance that describes the evolution of root-zone moisture.It represents the soil moisture as a continuity relationship between water inputs (snowmelt and rainfall minus interception), and outputs (evaporation and percolation to deeper layers) over several soil layers.The water balance is calculated separately for the three land surface types, each with a different number of layers.
Acknowledging that the evaporation of water from soil is mainly controlled by the available energy and the soil moisture conditions, final estimates of evaporation will be highly dependent on the reliability of the precipitation data driving the soil water budget.In order to constrain the resulting uncertainty in modelled evaporation, microwave remote sensing data of surface soil moisture are used.The running water balance estimates are corrected at the daily time step using a Kalman filter assimilation approach based on the estimated uncertainty of the satellite observations.

Inputs to the soil water budget
The inputs to the soil water budget come exclusively from precipitation, both as rainfall and as snowfall.Even though irrigation is not included as an input, the subsequent assimilation of the satellite soil moisture will partly account for it by adjusting the soil moisture seasonal dynamics of the area.
Precipitation is divided into rainfall and snowfall depending on the satellite observations of snow depth (D s ); when D s is over 10 mm (snow water equivalent), precipitation is considered snowfall (P s ).Rainfall (P r ) enters the soil directly, except for the fraction intercepted by tall canopies and evaporated back into the atmosphere (I ).P s however, does not enter the soil directly but accumulates in a layer on top of the soil column.This snow can either evaporate as E s (see Sect. 2.4), or melt and enter the soil water balance.The initial estimate of the snow depth (D − m ) for a given day (i) is calculated as This initial estimate is compared with D s .In the cases when the estimate exceeds the observed value, the difference is attributed to snow melt (F s ): and the estimated snow depth is reduced to match the satellite observation.The total flux of water into the soil water balance for day i is then calculated as In this study, the entire water flux (F ) infiltrates the soil column.With the intention of maintaining the simplicity of GLEAM, processes like surface overland flow (when the water flux exceeds the infiltration capacity of the soil) and bypass flow (when the flux reaches the aquifer directly) are considered to have a negligible effect in the evaporation processes at the coarse resolution of the model.Therefore, no horizontal movement of water or routing between adjacent pixels is considered in the methodology.

Root-zone water balance
In nature, the depth of the soil column that affects the evaporation rate depends on the rooting depth of the vegetation, and may vary from a few centimetres for grasses to as deep as four metres for forests.For bare soil, the lack of roots limits the thickness of the layer that affects evaporation to only a few centimetres.Because of those differences, the soil water balance is calculated for each land cover type individually.The shallowest soil layer has a depth of 0-0.05 m.For bare soil this is the only layer considered.For short vegetation a second layer is defined from 0.05-1.00m.Finally, for the fraction of tall canopy two extra layers are considered (0.05-1.00 m and 1.00-2.50m).
At each layer (l), the soil moisture content (w) on a given day (i) is modelled as: ig. 2. Schematic overview of the running water balance for the fraction of tall canopy (three layer rofile).In this example, the second layer is the wettest layer and therefore it determines the stress ctor, S (see Sect. 2.3).
37 Fig. 2. Schematic overview of the running water balance for the fraction of tall canopy (three layer profile).In this example, the second layer is the wettest layer and therefore it determines the stress factor, S (see Sect. 2.3).
where F (l−1) denotes the downward flux from the above layer, which in the case of the first layer will be the infiltration flux (F ) calculated through Eq. ( 4).E (l) represents the removal of soil water due to evaporation, z (l) is the thickness of the layer and F (l) is the percolation flux to the next layer.F (l) is estimated as the volume of water exceeding the field capacity (w fc ), hence The water percolating out of the deepest root-zone layer is assumed to be no longer available for plant uptake and therefore does not affect the modelled evaporation.Figure 2 presents an overview of the complete running water balance.

Satellite surface soil moisture assimilation
The depth of soil that affects microwave soil moisture observations is a direct function of the soil moisture conditions (see Ulaby et al., 1982).However, numerous studies in the past have shown that satellite-derived surface soil moisture is strongly related to the 0-0.05 m soil layer (e.g.Wagner et al., 2007;De Jeu et al., 2008;Draper et al., 2009;Gruhier et al., 2010).Therefore in GLEAM, satellite observations of soil moisture (θ) are assimilated with the modelled water content of the first soil layer (w (1) ) that is predicted by Eq. ( 5).The approach follows a one-dimensional Kalman filter design (see Crow, 2007).Every year, the methodology is first run without any data assimilation of soil moisture.Subsequently, time series of satellite observations at every pixel, are normalised to match the mean and variance of the time series of the model estimates with no assimilation.In addition, the cumulative density function of normalised satellite observations is scaled to match the cumulative density function of model estimates with no assimilation.The methodology is then run with assimilation of the scaled satellite observations.
The update of the model estimates at daily time step follows in which superscripts "−" and "+" denote values before and after the Kalman filter update.Krepresents the Kalman gain, which is calculated as where V denotes the error variance associated with the satellite observations (θ) and − is the background error variance of the Kalman filter forecasts.− is estimated as in which Q is the variance associated to the soil water balance estimates when propagated from time i − 1 to i. Then − is also updated as to obtain + , the variance error of the final estimates of soil moisture (w (1)+ ).
In our approach we consider a constant value of Q = 0.01.This implies that the value of K will be fully determined by the estimation of the variance error in the microwave observations (V ).According to De Jeu et al. (2008), the vegetation optical depth (τ ) can be used to approximate the polynomial relation existing between the uncertainty of the microwave soil moisture retrieval and the vegetation density.This relation can be described as The microwave soil moisture observations are obtained nearly every day when the temperatures are above freezing.Pixels covered by snow, presenting a fraction of open water larger than 20%, or those which show an annual negative correlation coefficient between time series of satellite observations and model estimates (with no Kalman filter) are not subject to this assimilation.In addition, satellite soil moisture is not assimilated in pixels presenting a fraction of tall canopy larger than 70%.This requirement is somehow redundant given the high value of τ in those pixels.The impact of this assimilation is explored in Sect.4.1 by comparison to in situ measurements of soil moisture.

Evaporative stress
For most of the land surface, the actual evaporation rarely -if at all -reaches the potential rate due to suboptimal environmental conditions.In those cases the actual evaporation will be less than the maximum rate for a given ecosystem.Environmental factors limiting the potential evaporation can be: a lack of available soil water, seasonal or occasional decrease  in biomass content, and extreme temperatures.To account for these effects it is common to define an empirical parameter (see for instance Barton, 1979) referred as evaporation stress factor (S), with unity indicating no stress, and zero indicating maximum stress.
In GLEAM, S is parameterised separately for tall canopies, short vegetation, and bare soil.This parameterisation is based on the soil moisture conditions, and (for the short vegetation fraction) a parameter accounting for the development of vegetation over the year (vegetation optical depth, τ ).
The soil moisture component of S is determined by the water content of the wettest soil layer as determined by the soil water module (see Sect. 2.2).This concept reflects the ability of vegetation to draw water from any layer within the root-zone, and affects the tall canopy (with three layers of soil) and short vegetation fractions (with two soil layers), but not the bare soil fraction (which presents only one layer of soil).For soil moisture values below wilting point (w wp ), the stress is the maximum (S = 0); for values above the critical moisture level (w c ), there is no stress (S = 1).Between w wp and w c the stress increases as soil moisture decreases following a parabolic function for the fraction of tall canopy, and an exponential relation for the fraction of short vegetation and bare land cover (Gouweleeuw, 2000;Owe and van de Griend, 1990).The stress functions used in GLEAM for the three land-surface components according to these parameterisations are defined and illustrated in Fig. 3.
The development of vegetation over the growing season as affected by environmental conditions and plant health is not modelled explicitly.Instead the microwave vegetation optical depth (τ ), is used as a proxy for the vegetation density because of its close relation to vegetation water content (De Jeu et al., 2008).In this study τ is used in short-vegetated covers to account for the effect of seasonal or occasional changes in biomass content on evaporation (i.e. because of harvesting, fires, etc.).Therefore, an important implication of using this dynamic estimate of vegetation density is that it adds variation to the otherwise static maps of cover fractions.
As an extra limit to the evaporative flux, the modelled evaporation is compared with the available water above w wp according to the soil water module (see Sect. 2.2).This assures no evaporation is extracted below w wp or from outside the root-zone.Priestley and Taylor (1972) showed that the Bowen ratio (the ratio of sensible to latent heat flux) would approach a constant value when air moves over a moist surface and gradients of temperature and specific humidity with height are small, or when the air becomes saturated with respect to moisture.The Priestley-Taylor (PT) equation has been shown to work well over many vegetation types with only small modifications.The formula calculates evaporation as a function of the available energy -net radiation (R n ) minus ground heat flux (G) -and a dimensionless coefficient (α) that parameterises the resistance to evaporation.Considering values of α for optimal environmental conditions (no evaporative stress), the model can be applied to describe the potential latent heat flux, λE p (MJ m −2 ), as:

Actual evaporation
where is the slope of the temperature/saturated vapour pressure curve (kPa K −1 ) and γ is the psychrometric constant (kPa K −1 ).λE p can be divided by the latent heat of vaporization, λ (MJ kg −1 ) -calculated as a function of temperature (Henderson-Sellers, 1984) -to derive potential evaporation (E p ) in mm.The magnitude of G is approximated in GLEAM as a fraction of R n , being 5%, 20% and 25% for the fraction of tall canopy, short vegetation and bare soil respectively (see e.g.Kustas and Daughtry, 2005;Santanello and Friedl, 2003).
For optimal environmental conditions (when actual equals potential evaporation), the value of α = 1.26 is welldocumented in the literature for grasslands.Similar values have also been found in past studies over bare land (Owe and Van de Griend, 1990;Caylor et al., 2005).However, Shuttleworth and Calder (1979) found that a value of α = 0.72 better reflected the conservative transpiration from forests; this As a result of suboptimal environmental conditions (due to soil water deficit or biomass changes), the volume of actual evaporation (E) is generally lower than the potential evaporation (E p ) calculated through Eq. ( 12).Several studies in the past (see for instance Barton, 1979) introduce the evaporation stress factor (S) to adapt the PT equation and account for the effect on E of suboptimal environmental conditions (see Sect. 2.3 for the parameterisation of S in GLEAM).In addition, when the canopy is wet the evaporation from forest is not well described by the PT equation (Stewart, 1977;Shuttleworth and Calder, 1979).In GLEAM, canopy rainfall interception is calculated independently (see Sect. 2.1).As a consequence of this separate estimation, the transpiration as calculated by Eq. ( 12) needs to be corrected by a fraction (β) of the interception loss (I ) to avoid the double counting of evaporation for those hours with wet canopy.Taking this correction into consideration, and adding the evaporation from the wet forest canopy and the effect of the evaporative stress, GLEAM describes E (in mm day −1 ) as: where β is considered a constant (β = 0.07 - Gash and Stewart, 1977).As mentioned before, E is calculated separately for the three land cover types, and subsequently aggregated to pixel scale through Eq. (1).For the fractions of short vegetation and bare soil, the I term in Eq. ( 13) is zero.Finally, the evaporation from snow-covered pixels is calculated by adapting and γ in the PT equation according to Murphy and Koop (2005).Literature values of α for snowcovered surfaces were not found and, therefore, α was calibrated based on 12 selected FLUXNET sites, each with more than fifty days of snow cover.It was found that α = 0.95 minimized the average error in cumulative sublimation for all sites.Moreover, snow-covered ecosystems are assumed to be unstressed due to the sufficient availability of water.Therefore in GLEAM, values of α = 0.95 and S = 1 are used as constants for every snow-covered pixel.

Satellite observations
The data used to run GLEAM in this exercise are listed in Table 1.All these data sets are primarily based on satellite observations.They are acquired from various sources and comprise well-validated products.Only the microwave vegetation optical depth represents a research product with limited validation.Its use in GLEAM for the parameterisation of the evaporative stress (Sect.2.3) and the estimation of the uncertainty of satellite soil moisture observations (Sect.2.2.3) is a unique feature of the proposed approach.The majority of the data sets are available at 0.25 degree regular grids; all the data sets presenting a different spatial resolution are re-gridded to a common 0.25 degree grid by means of Shepard's Method of inverse distance weighted interpolation (Shepard, 1968).

Net radiation
R n is the principal driver of the latent heat flux and the main input for the estimation of λE p by the PT equation (see Eq. 12).The NASA/GEWEX Surface Radiation Budget (SRB) Release-3.0contains global daily averages of surface longwave and shortwave radiative variables on a 1 • × 1 • grid.The data were obtained from the NASA Langley Research Center, Atmospheric Sciences Data Center NASA/GEWEX SRB Project.The product is based on a range of satellite instruments, reanalysis and assimilation.

Precipitation
The interception loss model (described in Sect.2.1) and the water balance (described in Sect.2.2) are driven by P as retrieved according to the Climate Prediction Center morphing Hydrol.Earth Syst.Sci., 15,[453][454][455][456][457][458][459][460][461][462][463][464][465][466][467][468][469]2011 www.hydrol-earth-syst-sci.net/15/453/2011/ technique (CMORPH) and provided by Joyce et al. (2004).This technique uses half-hourly infrared observations -Geostationary Operational Environmental Satellite (GOES), the Geostationary Meteorological Satellite (GMS) and Meteosat -to propagate higher quality microwave precipitation estimates from the Advanced Microwave Sounding Unit-B (AMSU-B), the Special Sensor Microwave Imager (SSM/I), the TRMM Microwave Imager (TMI) and the Advanced Microwave Scanning Radiometer (AMSR).Between measurements, intensity and shape of the microwaveobserved precipitation are modified by a time-weighted interpolation (morphing) resulting in a high spatial (0.07 • ) and temporal (30 min) resolution.Validation studies show better correlation between CMORPH and ground measurements than for most of the currently available satellite-derived precipitation products (Ebert et al., 2007).However, the spatial coverage of CMORPH is from 60 • N to 60 • S. In addition, it presents another practical disadvantage for its application in GLEAM: the product tends to underestimate precipitation at high latitudes, especially in winter-time (Zeweldi and Gebremichael, 2009;Tian et al., 2007).Consequently, for latitudes outside the CMORPH domain and snow-covered pixels, the 1 • daily Global Precipitation Climatology Project precipitation product (GPCP-1DD -see Huffman et al., 2001) is used.This product is produced by merging precipitation estimates from microwave, infrared, and sounder data observed by the international constellation of precipitation-related satellites, and precipitation gauge analyses (Huffman et al., 1997).GPCP-1DD has been widely used in different studies during the last few years as it represents one of the best available global precipitation products (Crow, 2007).
It is important to notice that the uncertainty in global precipitation products can be large according to the level of disagreement between the existing precipitation datasets.However, this uncertainty is difficult to estimate through comparison with gauge data due to the point-nature of precipitation ground measurements.In areas with dense observationalnetworks gauge-corrected products like GPCP-1DD are likely to outperform fully satellite-based products like CMORPH.The choice of CMORPH over GPCP-1DD for the exercise presented here, has to do with its high spatial resolution, full use of high quality TRMM observations and ability to capture orographic rainfall (see Hirpa et al., 2009).Neither CMORPH nor GPCP distinguish between rain and snow, and for this reason the satellite-observed snow depth (defined in Sect.3.3) is used to categorise precipitation as snowfall instead of the default classification as rainfall (see Sect. 2.2.1).

Microwave retrievals
An increasing number of geophysical land surface variables are successfully being retrieved from satellites carrying passive microwave radiometers.In general, microwave re-trievals have the benefit of being insensitive to clouds, resulting in a reliable twice-daily sampling rate.The methodology as presented here, relies on four of those variables derived from the AMSR-E radiometer on the AQUA satellite: surface soil moisture (θ), land surface temperature (T ), vegetation optical depth (τ ) and snow depth (D s ).The mean spatial resolution of the AMSR-E radiometer is between 12 km for the 36.5 GHz channel and 56 km for the 6.9 GHz channel.The first three parameters are derived with the Land Parameter Retrieval Model (LPRM) (Owe et al., 2008).LPRM is an iterative optimization and polarization index-based retrieval model that uses the dual polarization channels at a single low microwave frequency to derive θ and τ .In this paper the combined version (v04d) is used, in which the default 6.9 GHz based retrieval is replaced by the 10.7 GHz based product in areas that suffer from high levels of radio frequency interference in the lower band.
The LPRM soil moisture product has been validated over different ecosystems and is estimated to have an average accuracy of 0.06 m 3 m −3 (see De Jeu et al., 2008).Daily maps of satellite soil moisture used in GLEAM are derived from the next day's descending overpass (0130 LST) of AMSR-E.No gap-filling is applied.
The microwave vegetation optical depth presents a direct relation with vegetation water content (Kirdiashev et al., 1979).Here, a five-day central moving average is calculated to gap-fill the LPRM-derived τ .All LPRM products are limited to the non-frozen land surface, and therefore τ can still present long gaps in wintertime despite the five-day central moving average.These long gaps are filled with the 10th percentile of the values measured in a specific grid cell over the year.As τ represents a pixel-averaged value, it needs to be reassigned to each of the three land cover fractions.In GLEAM this is done based on two assumptions: (1) τ for the bare soil fraction is zero, and (2) τ for the short vegetation fraction is 60% of that of the tall canopy fraction.This 60% is based on the observed differences between values over entirely forested and nearby short-vegetated pixels.
LPRM uses the Ka-band (37 GHz) vertical polarised channel to retrieve the physical temperature of the emitting surface (skin temperature), a method recently described by Holmes et al. (2009).For vegetated surfaces this will be the temperature of the top of the canopy, a close estimate of the temperature required in the PT equation (see Priestley and Taylor, 1972).Daily maps of temperature used in GLEAM are derived as an average of the descending (01:30 a.m.) and ascending (01:30 p.m.) AMSR-E overpasses.A five-day central moving average is applied to gap-fill the data.Under frozen conditions the temperature is not retrieved from microwave data so gaps can still occur after the five-day central moving average.These long gaps are filled with air temperature data from the International Satellite Cloud Climatology Project (ISCCP) (Zhang et al., 2004).
Finally, the strong effect that snow has on the microwave emission is used by the National Snow and Ice Data Center (NSIDC) to retrieve snow depth.In this study we use the AMSR-E/Aqua daily L3 global snow water equivalent EASE-Grids V001 (Kelly et al., 2003).

Static data sets
A limited number of static data sets are used in the methodology.The most important one is the global Vegetation Continuous Fields product from MODIS, MOD44B (Hansen et al., 2005) which describes every pixel as a combination of its fractions of tall canopy, short vegetation and bare soil.The global fields of wilting point, critical soil moisture and field capacity are derived from the Global Gridded Surfaces of Selected Soil Characteristics (IGBP-DIS) (Global Soil Data Task Group, 2000).For the interception loss model, information to determine the mean rainfall rate is derived from the Combined Global Lightning Flash Rate Density monthly climatology from NASA (Mach et al., 2007).Finally, a digital elevation model is used to calculate the air pressure as it varies with height above sea level according to the barometric formula (and in accordance with the standard atmosphere).

Validation and discussion
The results presented here correspond to the application of GLEAM for the year 2005.A two-year period (2003)(2004) is used to spin up the soil water module.Both the soil moisture profile and the final estimates of evaporation are validated using in situ measurements.This exercise is complementary to the independent validation of the GLEAM interception loss estimates presented by Miralles et al. (2010).

Soil moisture profile validation
In situ measurements of soil water content from a selection of stations from the Soil Climate Analysis Network (SCAN -see Schaefer et al., 2009) are used to validate the daily soil moisture profile as modelled for the corresponding pixels.SCAN stations present soil moisture sensors at depths of 0.05, 0.1, 0.2, 0.5 and 1 m.Only SCAN stations with continuous measurements during the year 2005 are selected for this validation.These stations are located in grasslands or other short vegetation ecosystems within the US.To be consistent with the land use at the stations, only the modelled soil moisture for the fraction of short vegetation within the corresponding pixel is used in this validation.The Pearson's correlation coefficients between daily-averaged in situ measurements and modelled soil moisture content for the root-zone layers 1 and 2 (w (1) and w (2) respectively) are calculated at each station.Estimates of w (1) are compared with ground measurements at 5 cm; w (2) is compared with the average of the measurements at 0.05, 0.1, 0.2, 0.5 and 1 m.
Table 2 describes the stations used in this study as well as the individual correlations found between in situ measurements and GLEAM estimates of soil moisture.The mean correlation coefficients for a total sample of 30 stations are 0.60 and 0.69 for the first and second layer respectively.Soil moisture in shallow layers presents faster dynamics and is therefore less dependent on long term seasonal variations; this explains the slightly lower correlations found for the first layer of soil.The histogram for the first layer is presented in Fig. 4a, which also illustrates the effect of the assimilation of θ into the profile.Figure 4b shows the same inferences but for the second layer of soil.
An increase in the correlation in the first layer of soil is shown for 21 of the 30 stations when θ is assimilated.For this layer, the average correlation coefficient increases from 0.57 to 0.60.Even though the second layer is not subjected to the assimilation scheme, the Kalman filter update of w ( 1) is likely to have an impact on today's S.This may affect tomorrow's root-zone moisture profile not only by altering the initial w (1) but also by changing the volume of water removed from the profile through E (see Sect. 2.4).However, the improved characterization of w ( 1) is shown to have little effect on the time series of w (2) .This is mainly related to the fact that the much lower thickness of the first layer makes variations in this layer cause only subtle changes in the rest of the profile.In the near future, the assimilation of θ will also be done at deeper layers to propagate the effect of this assimilation more effectively.

Selection of ground stations
The modelled evaporation for the year 2005 has been compared with eddy covariance measurements at a sample of FLUXNET stations.FLUXNET is a global network of micrometeorological towers (see Baldocchi et al., 2001) with the principal aim of quantifying carbon, water vapour and energy fluxes.At each station the evaporation flux is measured using the eddy covariance technique, which samples a distance of 100 to 2000 m upwind of the tower.Given that the method is generally unreliable during rainfall, for this validation exercise we compare the modelled E without the I component (note that Miralles et al. (2010) already validated the GLEAM interception loss product against a set of independent mass balance evaporation measurements).
FLUXNET stations are mainly located in Europe and the US, but cover the most common vegetation types and climates.For the purpose of this validation a station by station quality check was performed based on: (a) the amount of gap-filling in each daily aggregate (only days in which less than 10% of the half hourly data to form the aggregate were gap-filled), (b) the subsequent availability of daily data for the study period (only stations with a coverage of at least 60% of the days in 2005), and (c) the quality of their energy balance closure (only stations with less than 50% mismatch in their energy closure).This yielded a total of 43 reliable FLUXNET sites covering a large variety of ecosystems.In the analysis below, these 43 stations are grouped based on the type of vegetation cover (short vegetation or tall canopy) and the volume of annual precipitation for the year 2005 according to CMORPH (dry: P <= 500 mm, wet: P > 500 mm), resulting in four functional groups.Therefore we distinguish between group: (A) tall canopy and wet climate (N = 10 stations), (B) tall canopy and dry climate (N = 9), (C) short vegetation and wet climate (N = 13), and (D) short vegetation and dry climate (N = 11).Table 3 presents the list of the 43 stations and their corresponding groups for the validation exercise.

Point versus pixel aspects
The FLUXNET observations are essentially point measurements when compared to the 0.25 degree resolution pixels of GLEAM-modelled E. As the methodology considers dif-ferent surface types (short vegetation, tall canopy and bare soil), it accounts for sub-pixel heterogeneity to a certain extent.In this validation analysis the ground observations are compared with the modelled E estimates corresponding to the specific land-surface type associated with the site.Despite this sub-pixel heterogeneity, it is important to notice that the input data of GLEAM consist of uniform values for the whole grid box.This is crucial when it comes to R n .The resolution of R n is the lowest (1 degree) of all primary input data, and the energy budget is highly dependent on the particular characteristics of the surface (e.g.albedo) and therefore on the land cover.Moreover, the weight of R n in the PT equation guaranties the propagation of these uncertainties and makes R n the most important input in the estimation of E p through GLEAM (in wet areas presenting values of S close to 1, R n will be responsible for the majority of uncertainty in the final E estimates).Therefore, in order to better compare the relative merits of the evaporation methodology over different vegetation types -and reduce the magnitude of the uncertainties related to the driving data -we also report the results of a model run that substitutes the stationmeasured R n for the satellite-based R n .

Time series validation
The statistics of the validation of the daily time series of E are summarized in a Taylor diagram (Taylor, 2001) in Fig. 5a.
For each of the four groups described in Sect.4.2.1, this figure displays the average correlation coefficient, standard deviation and RMSD of the comparison with the stations of the group.Both standard deviation and RMSD are normalised using the corresponding station as a reference, and therefore the point denoted as "Ref" represents the location in the diagram of the time series of every station.The origin of the arrows indicates the results using the satellite-based R n as input and the point of the arrows indicates the statistics with the site-measured R n as input.The use of in situ R n leads to a general improvement in the correlation and a reduction of the magnitude of the residuals.As expected, this improvement is unambiguous in wet regions (in which evaporation is mainly determined by the available energy), but the correlation coefficients increase for the four groups.In group A, the slight overestimation of the variance is also corrected; in group B however, the change has a negative impact in the variance of model estimates.For the remaining of the validation exercise, only the results of the run using the site-measured R n are considered.
In the second Taylor diagram (Fig. 5b) the results of the validation of E estimates at daily time step are compared with the results of the monthly aggregates.Unsurprisingly, the aggregation of the daily evaporation over the whole month results in an improvement of the model statistics, especially in terms of correlations.In group A this improvement is more subtle due to the small amplitude of the seasonal cycle found in tropical forests -the station in Amazonia is the only one of the 43 stations that shows degradation in R. As it can be appreciated in Table 3 (which presents the values of the correlation coefficients for the individual locations in the two right columns), the Amazonian site on its own is responsible for the lower average correlation coefficient for group A found in Fig. 5. Stations in group C present a high average correlation with FLUXNET data (R = 0.85 for the daily and R = 0.91 for the monthly time series) in agreement with the original intention of the Priestley-Taylor method to estimate evaporation from short unstressed vegetation.Despite the fact that groups B and D show the highest correlations, it is important to note that in dry regions the presumed larger amplitude of the seasonal cycle of evaporation is likely to have a positive effect on the correlation coefficients.Overall, there is a high correspondence of GLEAM estimates with FLUXNET observations for each of the four groups, both at daily and monthly time-step.The average correlation for the 43 stations is R = 0.83 for the daily and R = 0.90 for the monthly series.Both the transpiration from tall canopies and short vegetated ecosystems seem to be equally well characterised by the model.Moreover, the extra complexity introduced by the modelling of evaporation stress does not seem to have a negative effect in the performance of the methodology over dry regions.Figure 6 presents an example of the daily time series of GLEAM estimates and FLUXNET observations for each of the four groups.The Amazonian station was selected to show how the lack of a clear seasonal cycle can affect the statistics in the comparison (and therefore the results for group A in Fig. 5).The other three stations are representative examples of each of the three groups they belong to, and show the good correspondence between in situ and model estimates of evaporation over different ecosystems.
In the last few years, FLUXNET data have been used to evaluate other methodologies dedicated to derive global evaporation estimates from satellite observations.Zhang et al. (2010) did a similar comparison to the one presented here.They also reported a good correspondence between model estimates and in situ observations.However they did not present the average value of the correlation coefficients at the stations, but the correlation coefficient that resulted from plotting the estimates from all the stations together in one single scatter-plot.Fisher et al. (2008), on the other hand, listed the individual correlations found at every station in their comparison of monthly estimates and station-based monthly aggregates of evaporation.The results of the validation of the methodology proposed by Fisher et al. (R = 0.93, N = 16) are comparable to the ones found here for the validation of GLEAM (R = 0.90, N = 43).

Annual totals and bias
With the aim of identifying a possible systematic bias for any of the four groups, Fig. 7 compares the total volumes of modelled and measured E for 2005 at each of the 43 FLUXNET stations (these volumes are also presented in Table 3).The correlation coefficient shows a value of R = 0.80.The bias is as low as −5%, which represents an average underestimation of less than 20 mm yr −1 .These inferences only indicate the level of agreement between observed and modelled annual aggregates, and therefore they only show the skill of the model to capture the global distribution of annual evaporation (see next section).
Overall, none of the four groups presents a major bias; this indicates that the scatter seen in Fig. 7 is not likely to be a response to systematic errors in the parameterisation of the two different vegetation types, or the two climate conditions considered to define the groups.Nevertheless, the cumulative error at some of the stations can become important and it ranges between −48% to +73%.The standard deviation of the residuals is therefore high, as can be seen from the value of RMSE = 110 mm yr −1 .However, attending to the level of mismatch of the energy closure observed at some stations (see Table 3), FLUXNET evaporation measurements may also be greatly biased.Consequently a validation exercise based on the correlations at each station individually (like the one performed in Sect.4.2.3) may be a better indicator of the performance of the methodology.

Global application of the methodology
The map of evaporation for 2005 as modelled by GLEAM is presented in Fig. 8.The spatial patterns appear reasonable and the range of values corresponds well with previous attempts to estimate global evaporation (see Jiménez et al., 2011).A detailed study of the spatial distribution of the GLEAM-modelled E is the topic of a further paper that will analyse the global magnitude of the latent heat flux and its seasonal variability, the relative importance of rainfall interception loss, the global distribution of water available for runoff and the physical processes controlling transpiration over the different regions of the world.

Conclusions
Evaporation remains one of the biggest unknowns within the global water balance.Improved representation of its global dynamics is essential to lead to a better understanding of the expected acceleration of the hydrological cycle.There have been several recent efforts towards the development of observation-based estimates of global evaporation; these attempt to create independent, daily-data driven benchmarks for GCM developers to improve their predictions of future climate.GLEAM (Global Land surface Evaporation: the Amsterdam Methodology) represents a new approach that combines Hydrol.Earth Syst.Sci., 15, 453-469, 2011 www.hydrol-earth-syst-sci.net/15/453/2011/ a wide range of currently existing satellite-sensor products to estimate reliable fields of daily global evaporation at a 0.25 degree spatial resolution.Because the methodology is based on the Priestley and Taylor (1972) radiation-driven evaporation model, it limits the number of spatially-varying surface fields that need to be specified and cannot be detected from space.The applicability of GLEAM relies exclusively on the availability of a suite of remotely-sensed input data products.Its simple strategy allows the application of the methodology, not only at a global scale (i.e.studies of trends in evaporation, evaluation of GCMs' performance, etc.), but also at a watershed scale through the utilisation of better resolution input data (i.e.radiometers, in situ observations, etc.).Its minimal dependence on static fields of variables -unlike many other models -avoids the need for parameter tuning and makes the quality of the evaporation estimates rely on the accuracy of the satellite inputs.As satellite-based observations are not error-free, the approach could potentially benefit from the assimilation of in situ observations in areas with dense ground-observational networks.A major distinguishing feature of the methodology is the detailed estimation of satellite-derived global fields of forest rainfall interception.Other characteristics are the coupling of the radiation-driven transpiration to the ground bio-physical processes (due to the parameterisation of the root-zone evaporative stress condition), and the separate estimation of bare soil evaporation and snow sublimation.
Model estimates have been successfully compared with ground data from a wide range of ecosystems.The two main intermediate products of GLEAM have been individually validated: the forest rainfall interception (R = 0.86, Bias = −0.6%,N = 42 -in Miralles et al., 2010) and the root-zone soil moisture (R = 0.60 and R = 0.69 for surface and deep layers respectively).In addition, final evaporation estimates have been validated against one year of eddy covariance measurements from 43 FLUXNET stations.Results show a high average correlation with ground measurements, both at a daily (R = 0.83) and a monthly (R = 0.90) time scale.Moreover, no systematic bias for specific vegetation types or rainfall conditions has been detected.
Updates to the methodology are planned in the assimilation of remotely-sensed soil moisture data.These updates include the characterisation of the variance of soil water balance estimates (Q), and the assimilation of satellite observations into deeper layers to better propagate the optimisation through the entire root-zone.
In an ongoing study we analyse the spatial distribution and magnitude of the global estimates of latent heat flux, and their seasonal variability and relative importance of their   different components; this includes an insight into the global distribution of the evaporation drivers and the generation of water available for runoff.Our ultimate goal is to extend the time period to produce a global 0.25 degree daily evaporation data set spanning from 1984 to present.Considering the availability of the different input data sets over time (presented in Table 1), this exercise will require the use of different precipitation and snow water equivalents products.The extended evaporation data set will be compared with other existing products in forthcoming studies integrated within the LandFlux-Eval initiative (Jiménez et al.,2011;Mueller et al., 2011).GLEAM products will be made available in the VU University Amsterdam geoservices website: http://geoservices.falw.vu.nl.

Fig. 3 .
Fig. 3. Overview of stress parameterisations for tall canopy, short vegetation (at two levels of vegetation density: (a) τ = 0.2, (b) τ = 0.8), and bare land.The values of w wp and w c are considered to be 0.1 and 0.3 m 3 m −3 respectively; w w corresponds to the soil moisture modelled for the wettest layer.

Fig. 4 .Fig. 4 .
Fig. 4. Histograms of the correlation coefficient (R) of the modelled soil water content with in situ SCAN data for: (a) first layer of soil, (b) second layer of soil.The histograms show the difference in the validation with and without the data assimilation (DA) of satellite soil moisture.

Fig. 5 .Fig. 5 .
Fig. 5. Taylor diagrams of the validation results for the groups listed in Table 3. RMSD and standard deviation are normalised against the reference represented by the time series of the corresponding FLUXNET station; therefore, the point denoted as "Ref" represents the location in the diagram of the time series of every station.(a) shows the results of the comparison between daily time series of modelled E and the E measured at the 43 FLUXNET stations.The dots correspond to the statistics of the model run with the satellite R n as input; the arrows point the results of the model validation when using the R n measured at the stations as input.(b) shows how the statistics improve when comparing monthly averages instead of daily time series (using the station R n as input).

Fig. 6 .
Fig. 6.Examples of daily time-series of FLUXNET and GLEAM E (in mm) from each of the four groups defined in Sect.4.2.1

Fig. 7 .Fig. 7 .
Fig. 7. GLEAM annual E against annual cumulative evaporation at the 43 FLUXNET stations for th year 2005.Stations are grouped by vegetation cover and climate conditions (see Sect. 4.2.1).

Table 1 .
Remotely-sensed data sets used for computing GLEAM E estimates (see Sect. 3 for explanation of abbreviations).Satellite value was estimated for two forest stands in the UK, where soil moisture deficit could be considered low although no parameterisation of the stress due to soil moisture conditions was performed.In 1984, Shuttleworth et al. found that a value of α = 0.91 better suited the parameterisation of forest potential evaporation in a tropical region.In GLEAM, a constant value of α = 0.8 is used to parameterise the tall canopy fraction, while a value of α = 1.26 is applied in both the short vegetation and bare soil fractions.

Table 2 .
SCAN study sites and results of the validation of the modelled soil moisture profile.

Table 3 .
FLUXNET sites used in the validation.The coefficient e represents the mismatch in the energy balance at the station.The measured and modelled cumulative evaporation (for 2005) and the correlation coefficients of the comparison at daily and monthly time-step are also listed (in situ R n used).