Hydrology and Earth System Sciences A Bayesian approach to estimate sensible and latent heat over vegetated land surface

Sensible and latent heat fluxes are often calculated from bulk transfer equations combined with the energy balance. For spatial estimates of these fluxes, a combination of remotely sensed and standard meteorological data from weather stations is used. The success of this approach depends on the accuracy of the input data and on the accuracy of two variables in particular: aerodynamic and surface conductance. This paper presents a Bayesian approach to improve estimates of sensible and latent heat fluxes by using a priori estimates of aerodynamic and surface conductance alongside remote measurements of surface temperature. The method is validated for time series of half-hourly measurements in a fully grown maize field, a vineyard and a forest. It is shown that the Bayesian approach yields more accurate estimates of sensible and latent heat flux than traditional methods.


Introduction
Sensible and latent heat (i.e.evapotranspiration) fluxes between the land surface and the atmosphere are important components of the land surface energy balance.Different techniques exist to estimate them, generally based on micrometeorological methods, including the eddy-covariance, Bowen ratio technique and bulk transfer equations, for example.Alternatively, estimates of evapotranspiration can be obtained from the soil water balance (see e.g.Verhoef and Correspondence to: C. van der Tol (tol@itc.nl)Campbell, 2005, for an overview of both types of methods), from which sensible heat flux could be derived, if values of net radiation and soil heat flux were available.
Field scale measurements are used to study the local water and energy balance and to gain process understanding, but they are often not representative of large areas.Remote sensing measurements provide a spatial coverage, but not all variables needed to estimate sensible and latent heat can be measured by remote sensing.
It is common practice to combine remote sensing and field data to estimate sensible heat, H , and latent heat flux, λE, spatially (e.g.Su, 2002;Kustas et al., 2007;Anderson et al., 2008).Remote sensing products used as input include brightness temperature and emissivity (Bastiaanssen et al., 1998), reflected shortwave radiation and NDVI, which are used for the derivation of vegetation structure and aerodynamic resistance (Su, 2002).Weather station variables include wind speed, air temperature and humidity.Latent heat flux is either solved as a residual term in the energy balance, with H obtained from a bulk transfer equation or directly calculated using the bulk transfer equation for latent heat flux.The problem with the latter approach, however, is that estimates of specific humidity at the land surface are hard to obtain.
There are various sources of uncertainty associated with these flux estimates based on remote sensing and field data: the representativeness of the weather station data, the atmospheric correction of satellite data, the derivation of remote sensing products, the relationship between remote sensing products and surface characteristics (such as roughness length, which affects aerodynamic conductance), the model itself, and the interpolation in time between satellite overpasses.
Published by Copernicus Publications on behalf of the European Geosciences Union.Air and surface temperature T a , T s K Vapour concentration in the air, and in soil and leaf pores q a , q s kg m −3 Quantification of the uncertainty is an essential part of surface energy balance modelling.One of the techniques to handle uncertainty of data is data assimilation (Troch et al., 2003).Data assimilation makes it possible to quantify uncertainty, and to find the best solution of a model given various sources of input data with associated errors.Time series used for surface energy balance modelling are often incomplete, and data are collected at different spatial scales.In such cases, data assimilation is a useful technique.It has been used to improve the predictions of hydrological models (Schuurmans et al., 2003), to retrieve spatial and temporal trends in root water extraction by vegetation from remote observations of surface temperature and ancillary data (Crow and Wood, 2003), and to retrieve (simulated) soil moisture by model inversion (Reichle et al., 2002).
This paper presents a data assimilation approach to estimate sensible and latent heat flux over vegetated land surface.The data used to test the approach include remotely measured surface temperature in conjunction with ground-based meteorological measurements.In principle, the meteorological data are sufficient to calculate sensible and latent heat flux.This makes it possible to calculate the energy fluxes even in the absence of surface temperature measurements, for example between satellite overpasses or during cloudy conditions.By adding remote measurements of surface temperature (by means of data assimilation), the estimates of the fluxes are improved, because information is added to the model.This approach has already been suggested by Franks et al. (1999).A simple form of data assimilation is used in this paper: the Bayesian theorem described in statistical handbooks (e.g.Carlin and Louis, 1996).
In Sect.2.1 we first present the energy balance model, bulk transfer equations, and discuss classic approaches to estimate evapotranspiraton (Sect.2.2).Next, the Bayesian theorem is explained in detail (Sect.2.3).Furthermore in Sect. 3 it is described how the data are used to test the Bayesian approach, and the field experiments are described.We finish with the results and discussion in Sects.4 and 5, respectively.

Energy balance
The standard energy balance equation is used, consisting of four components: net radiation, R n , soil heat flux, G, sensible heat flux, H and latent heat flux, λE (all in W m −2 ).Energy involved in melting and freezing or in chemical reactions, energy stored in the canopy, and energy horizontally transported by advection are ignored.Hence (Brutsaert, 2005): Parameter λ is the latent heat of evaporation of water (J kg −1 ) and E the evapotranspiration rate (kg m −2 s −1 ).
Sensible and latent heat flux are the fluxes of heat and water vapour between the surface and the atmosphere carried by turbulent air flow.They are calculated as the product of a conductance and a driving force between the surface (subscript "s") and the air at measurement height (subscript "a"), often called "bulk-transfer equations": (2) where ρ is the mass density of air (kg m −3 ), c p the heat capacity of air (J kg −1 K −1 ), g a is the aerodynamic conductance (m s −1 ) for transport of heat and vapour from the surface boundary layer into the atmosphere, g s (m s −1 )is the surface conductance for water vapour transport from stomatal cavities in leaves or soil pores, to the leaf or soil surface boundary layer, T s is surface temperature and T a is air temperature (both in K), q s is the vapour concentration in stomata or soil pores (kg m −3 ), and q a is the vapour concentration in the air (kg m −3 ).
Equations ( 1) to (3) contain 13 parameters and variables.Three of them can be considered constant or can accurately be estimated (ρ, c p and λ).Out of the other 10 (Table 1), at least 7 need to be measured or estimated in order to solve the three equations.Weather station data are commonly used for the driving variables, such as T a , q a or wind speed (required for calculation of r a ); insofar they are representative for the area under study.
The surface conductance and the aerodynamic conductance are particularly difficult to estimate.
Surface conductance, g s , of individual leaves can be derived from leaf gas exchange measurements, but scaling to canopy level is not trivial (Baldocchi et al., 1991).Surface conductance varies with the actual photosynthesis rate, soil moisture content, air humidity, air temperature.At canopy level, empirical and semi-empirical relationships between surface conductance and environmental conditions are commonly used (Jarvis, 1976;Cowan, 1977;Ball et al., 1987;Hydrol. Earth Syst. Sci., 13, 749-758, 2009 www.hydrol-earth-syst-sci.net/13/749/2009/ Leuning, 1995).These relationships require vegetationspecific, a priori coefficients or local calibration against measured fluxes of water and carbon dioxide.However, for many applications these models are too detailed, as parameters for the vegetation under study are often not available.The aerodynamic conductance, g a , depends on the momentum and scalar roughness length of the surface, z 0M (m) and z 0H (m), respectively, wind speed u (m s −1 ), and stability of the atmosphere, as combined in the logarithmic wind profile (Tennekes, 1973;Garratt, 1992).The momentum roughness length, z 0M , is usually estimated from leaf area index and vegetation height (e.g.Raupach, 1994;Verhoef et al., 1997a), whereas the scalar roughness length, z 0H , is often taken as a constant fraction of z 0M ; this fraction varies widely, largely depending on canopy density (Stewart et al., 1994;Verhoef et al., 1997b).Optical remote sensing is often used to estimate vegetation height and leaf area index (Su, 2002).However, estimates of aerodynamic resistance based on optical remote sensing contain a large degree of uncertainty.
Apart from the aerodynamic and surface conductance, the difference between surface and air temperature (T s −T a ) is also difficult to estimate.Sensible heat flux is proportional to this difference; a relatively small difference between two measurements.When T s is taken from remote sensing and T a from weather station data, then errors in either T s or T a may cause significant errors in the estimates of H .

Classic approaches to estimate evapotranspiration
The difficulty in estimating g a , g s and (T s −T a ), in Eqs.(1) to (3) has been overcome in a number of ways.This has resulted in three techniques to calculate evapotranspiration from the energy balance: 1. FAO-approach: evapotranspiration is calculated by multiplying a "reference evapotranspiration" by an empirical coefficient for a specific crop, extracted from a table based on expert knowledge.The reference evapotranspiration is calculated for a crop with known values for g a and g s , usually the typical short, well-watered grass of a meteorological station.This method is disseminated by the FAO (Allen et al., 1998(Allen et al., , 2006)), and is mainly used to calculate irrigation requirements.
2. T s -approach: evapotranspiration is calculated by solving Eqs.
In that case, on top of the standard meteorological variables, surface temperature measurements and estimates of aerodynamic conductance g a are required.This technique is used in remote sensing, for example in the model SEBS (Su, 2002).Surface temperature is retrieved from thermal remote sensing, and aerodynamic resistance from optical remote sensing.
3. Conductance-approach: evapotranspiration is calculated by solving Eqs. ( 1) to ( 3) with H , λE and T s as unknowns.Aerodynamic conductance is estimated from a measured wind profile or from a priori values for a specific vegetation type, and surface conductance g s is estimated with one of the combined photosynthesisconductance or empirical models discussed above.Values for the parameters of such models are found from calibration of the model against measurements of g s or taken from the literature.

Bayesian approach
The Bayesian approach used in this paper, as an alternative to the classic approaches described in Sect.2.2, combines the methods ( 2) and (3) discussed above.In the Bayesian approach, neither a real distinction between "knowns" and "unknowns" is made, nor a distinction between "variables" and "parameters".The distinction is rather between "a priori" variables and parameters (values for parameters and variables that were measured, estimated, or taken from the literature), and "posterior" values (most likely values for parameters and variables, taking into account all a priori values and their uncertainty, and using the model).A priori values are indicated with a tilde (∼) placed above a symbol.Posterior values are indicated with a hat (ˆ).The "real" values, which are never known, do not have a superscript.
It is now assumed that besides ρ, c p and λ, also T a , q a , R n and G are accurately measured, such that their uncertainty can be ignored.For g a , g s and T s , uncertain, a priori estimates are used, and for q s , the saturated humidity at T s is used.No a priori information about the fluxes H and λE is used.Doing so results in three equations (Eqs. 1 to 3) with two unknowns (H and λE), which implies that the energy balance is over-specified.Theoretically, Eq. ( 1) is not no longer strictly necessary (we could estimate the fluxes without using Eq.1), but in doing so we abandon the physically realistic and useful constraint of energy balance closure.Similarly, we could calculate the fluxes without using measured T s (as in the conductance-approach), but in doing so we ignore part of the data that we have available.Finally, we could calculate the fluxes without using a priori g s (as in the T s -approach), but in doing so we ignore that g s is limited to physically realistic bounds.The fact that the model is over-specified is not a problem, since the values for g a , g s and T s are not fixed but uncertain.They will be adjusted, resulting in posterior values.
In the T s -approach of Sect.2.2, all uncertainty of the measurements propagates into the unknown g s , which may lead to unrealistic values for g s .Similarly, in the conductanceapproach of Sect.2.2, all uncertainty of the measurements propagates into the unknown T s , which may lead to unrealistic values for T s .In the Bayesian approach, all three parameters g a , g s and T s are constrained.In this way, unrealistic values for any of them are avoided.It is expected that this also leads to more accurate estimates of the fluxes H and λE.
The model (Eqs. 1 to 3) is re-written such that T s is output of the model, by eliminating H and λE from Eqs. (1) to ( 3) and writing T s as the dependent variable.In other words, T s is a function of all input variables and parameters.Hence: where x represents the measured output (in this case measured T s ), f the model equations, θ a vector of a priori values of variables and parameters, and w background noise caused by the uncertainty of the model and measurements.Note that Eq. ( 4) could be expressed in terms of the parameters of Eqs. ( 1) to (3) explicitly, but this would result in a rather long expression.The vector notation is used if the problem is solved for multiple time steps, multiple pixels, or for multiple measurements and parameters.This aspect is discussed later in this paper.
Parameter values (g a and g s ) have an a priori probability density function p(θ), and measurements of T s a probability density function p ( x θ).Note that g a and g s are variables.In this study we will calculate them for each time step independently.For an individual time step, g a and g s can be considered as model parameters.The probability density function of the parameter values is determined by the accuracy of a priori estimates, whereas the probability density function of the measurements of T s is determined by the accuracy and representativeness of the measurements.The implementation of these probability density functions for the data used in this study is discussed in Sect.3. The posterior probability density of the parameters (i.e. the parameter values, given the measurements of T s ), is calculated with the classic Bayes' theorem (Carlin and Louis, 1996): The numerator in Eq. ( 5) is a multiplication of the two probability density functions, whereas the denominator is a normalization term.A low probability of a value for T s (a T s that deviates significantly from the measured value), or a low probability of a value for the parameters (a conductance that deviates significantly from the a priori value), results in a low probability of the posterior parameter values.The expected values of the parameters are given by the minimum least square error estimate of the parameters.These expected values can be calculated by integrating the probability density function: The expected values are the posterior (or the optimum) values for the parameters.The posterior values are the weighted average (weighted over their probability) of all possible parameter values.The integration in Eq. ( 6) is avoided by calculating only the peak of the probability density function: where C x and C θ are the covariance matrices for the measurements and the parameters, respectively.Parameters θ are the posterior parameter values for g a and g s used later to calculate sensible and latent heat flux.
The last part of Eq. ( 7) (between the brackets) is a cost function or the quadratic error.The quadratic error is the sum of two components: the quadratic error in the parameters (θ) scaled with the uncertainty of the a priori estimates, and the quadratic error in the measurements of T s , scaled with the uncertainty of the measurements of T s .The optimum parameter values are located at the minimum total error.If both p( x |θ ) and p(θ) are Gaussian, then the solution is exact (and Eqs. ( 6) and ( 7) will give equal results), otherwise it is an approximation.In this study, Gaussian distributions for both functions are assumed.
The key issues are to estimate a priori values of θ, and to estimate the two covariance matrices.These matrices describe the uncertainty of all input (measurements and the a priori values), and determine the contribution of different input variables to the posterior estimates.These issues are addressed in detail in Sect.3.1.

Model input
The Bayesian approach is applied using meteorological time series for three land cover types: maize, a vineyard and a forest, measured during intensive field campaigns (Sect.3.2).Measured values of four variables are directly used as input for the model: R n , G, T a , and q a .It is assumed that these values were accurately measured and that the measurements were representative of the land cover (R n = Rn , etc.).This implies that for these four variables, a priori and posterior values are equal to each other.In contrast, uncertainty is attributed to the other variables, g a , g s and T s .Posterior values for g a , g s and T s were calculated with the Bayesian approach (Eq.7).The posterior values were used to calculate sensible and latent heat fluxes with Eqs.(1) to (3).Measured sensible and latent heat fluxes were used for validation only.
It is necessary to estimate the a priori values for g a and g s , as well as the uncertainty of a priori g a , g s and T s .Hence, the following equations for g a , g s and T s are introduced:   where w indicates noise, θ 1 is a dimensionless parameter which includes the effects of surface roughness and stability of the boundary layer, u is wind speed (m s −1 ), and parameter θ 2 is surface conductance (m s −1 ).The expression for g a is an additional model, that enables us to estimate aerodynamic conductance a priori (g a is a function of wind speed).The a priori estimate for parameter θ 1 for neutral conditions is calculated using the logarithmic wind profile following Tennekes (1973) and Goudriaan (1977): where d=0.67h, z 0M =0.13h, z 0H =0.1z 0M , h is the vegetation height, z the measurement height of wind speed (all in m), and κ (=0.4) is Von Kármán's constant.Equation ( 9) does not include a correction term for non-neutral atmospheric conditions.For the a priori estimate of parameter θ 2 , the FAO standard value for short, well watered grass of g s =0.0143 m s −1 is used for all three study sites.Values for Ts are computed from Stefan-Boltzman equation, using measured outgoing longwave radiation and an emissivity of 0.98.For the forest site, no reliable measurements of outgoing longwave radiation were available.For this site, surface temperature was measured with small Negative Temperature Coefficient thermistors (NTC) attached to needles (referred to as "contact temperatures" in this paper).
The uncertainty of the a priori estimates and measurements, which determine the matrices C x and C θ in Eq. ( 7), are estimated as follows.It is assumed that the probability of θ 1 , θ 2 and T s have normal distributions with a standard deviation of one quarter of the difference between their upper and lower limits found in the literature.The upper and lower limit of θ 1 for crops (in this case maize and vineyard) and for forest are based on minimum and maximum values for z 0M reported in a review paper of Garratt (1993).For these extreme values of z 0M , corresponding extreme crop heights (h=z 0M /0.13) and zero plane displacement height (d=2/3h) are calculated.The upper and lower limits of θ 1 are calculated from the corresponding values of z 0M , z, h and d with Eq. ( 9).An empirical equation for g s of Allen et al. (1998) is used to estimate the upper and lower limits of θ 2 : where LAI active is the leaf area index (m 2 leaf m −2 surface) that contributes to transpiration.For the upper and lower limits of θ 2 (i.e.g s ), the values corresponding to a LAI active of 3.5 and 0 are used, respectively.The upper and lower limits of T s are estimated by applying Stefan-Boltzmann equation for two extreme values of emissivity (0.90 and 0.99 for the vineyard and 0.95 and 0.99 for the fully grown maize).For the forest site, the standard deviation of the readings of nine thermistors is used as a proxy for the standard deviation of T s .Table 2 presents the a priori values and standard deviations σ θ 1 , σ θ 2 and σ T s derived in this way.
It is further assumed that the covariances (cov(θ 1 , θ 2 ), cov(θ 1 , T s ), cov(θ 2 , T s )) are zero, and that the errors in θ 1 , θ 2 and T s of consecutive time steps of a time series are uncorrelated.The latter assumption may not be realistic: estimates of roughness parameters and surface temperature measurements may be biased and errors in θ 1 and T s are therefore most likely to be similar for consecutive time steps.However, it appears that reasonable results can be obtained even when these effects are ignored.Assuming that all covariances are zero makes it possible to solve Eq. ( 7) for every time step individually.This is computationally more efficient than solving Eq. ( 7) for the whole time series at once, which requires manipulation of large matrices containing the parameters of all (half-hourly) time steps.
The posterior parameters for an individual time step are now calculated with Eq. ( 7), using: The posterior (optimum) parameter set, θ, is calculated with Eq. ( 7).The minimum of the square error (the term between the brackets of Eq. 7) is found with the Nelder-Mead method.These posterior values are used to calculate sensible and latent heat flux with Eqs.(1) to (3).For comparison, sensible and latent heat fluxes are also calculated using the latter two methods presented in Sect.2.2.That is, Method 2: the T s -approach, using Eqs.( 1) to ( 3) with H , λE and g s as unknowns, and using the a priori expression for g a (g a = θ1 u) and measured T s ; and Method 3: the conductance-approach, using Eqs.( 1) to (3) with H , λE and T s as unknowns, and using a priori values for g a and g s .

Experimental setup
The Bayesian approach is applied using micrometeorological time series obtained for a maize crop (Sonning, UK), a vineyard (Barrax, Spain) and a Spruce forest (Speulderbos, The Netherlands).For the maize crop, 9 days, for the vineyard, 6 days and for the forest, 3.5 days worth of half-hourly data are used.

The maize field
Meteorological input variables and fluxes were obtained over a fully grown maize field (row spacing was 0.75 m, the within-row spacing was about 0.12 m) located at the Crops Research Unit at Sonning farm (a research facility owned by the University of Reading, UK).It is located 4 km from Reading in Sonning (UK), at 51.48 • N, 0.90 • W, elevation 35 m above sea level (see Houldcroft, 2004); the soil type is loamy sand.Data from 7 to 16 September 2002 were used, when the leaf area index was 3.4 m 2 m −2 and the canopy height 1.9 m.Net radiation was measured using a CNR1 fourcomponent net radiometer (Kipp and Zonen, Delft, The Netherlands) at a mounting height of 2.5 m above the ground.Soil heat flux was calculated using a Fourier analysis on measured soil temperatures, combined with estimates of thermal diffusivity and heat capacity, the so-called analytical or exact method (see e.g.Verhoef, 2004).The soil temperatures were acquired with type-NT 10 k thermistors (RS, UK) that had been encapsulated with a stainless steel housing, accuracy of ±0.2 • C, installed at nominal depths of 2 and 5 cm.Soil heat flux at the surface (z=0), i.e.G, was calculated by using a negative z (i.e.−0.02 m) in the analytical equation of soil heat flux.Thermal diffusivity was calculated using the Arctangent method (Verhoef et al., 1996), using soil temperature signals at both depths, and heat capacity was calculated from the soil moisture content at 5 cm depth, measured using a Thetaprobe (Delta-T Devices).
Wind speed was measured using an AN1 cup anemometer (Delta-T Devices, UK), air temperature and humidity were measured with a RHT2 psychrometer, all at a height of 4 m above the surface.Surface temperature was estimated by inverting Stefan-Boltzmann's law, using outgoing longwave radiation measured with the CNR1 radiometer and an emissivity of 0.98.Although contact temperature measurements of the surface were available as well, these were not used in order to approximate remote sensing measurements as much as possible.
Sensible and latent heat fluxes were obtained using a combination of a Solent R3 sonic anemometer (Gill Instruments Ltd, Lymington, Hampshire, UK) and a differential closedpath infrared gas analyser (LI-7500, LICOR Inc., Lincoln, NE, USA), with appropriate correction procedures applied.

The vineyard
Meteorological input variables and fluxes were obtained over a vineyard (row spacing was 3.35 m, the within-row spacing was approximately 1.5 m, LAI was 0.52 m 2 m −2 , fractional vegetation cover 0.33 and vegetation height about 2 m) located at the Barrax agricultural test site in Spain (39.06 • N, 02.10 • W), where various crops were grown, some of them irrigated.Data were collected between 14 and 21 July 2004 during an intensive field campaign (SPARC).The experiment has been described in detail by Su et al. (2008).
Net radiation was measured using a CNR1 fourcomponent net radiometer at a mounting height of 4.8 m above the ground.Soil heat flux was measured with 3 Hukseflux HFP01 heat flux plates (Campbell Scientific Inc., USA) at 0.5 cm depth.Heat storage above the sensors was neglected.Air temperature and relative humidity were measured with a HMP45 sensor (Vaisala, Finland) at 4.78 m height, wind speed with a cup anemometer (Vector Instruments, Ltd., United Kingdom) at 4.88 m height.Surface temperature was estimated by inverting Stefan-Boltzmann's law, using outgoing longwave radiation measured with the CNR1 radiometer and an emissivity of 0.98.All data were collected at 1 min intervals, and 10-min averages were stored, and half hourly averages were used in this study.
Sensible and latent heat fluxes were obtained using a combination of a Solent R3 sonic anemometer (Gill Instruments Ltd, Lymington, Hampshire, UK) and a differential closedpath infrared gas analyser (LI-7500, LICOR Inc., Lincoln, NE, USA) installed at 3.4 m height.Further details about the experiment can be found in Timmermans et al. (2009).

The forest
Meteorological input variables and fluxes were obtained over a Douglas fir stand planted in 1962 in the Veluwe forest ridge in the Netherlands (52.50 • N, 05.69 • E).The research site is equipped with a 47 m high measurement tower maintained by the Dutch National Institute for Public Health and the Environment (RIVM).The tree density is 785 trees per hectare and the tree height 32 m.Leaf area index is approximately 5 m 2 m −2 .The topography is slightly undulating with height variations of 10 to 20 m within distances of 1 km.Data were collected between 10 and 21 June 2006 during an intensive field campaign (EAGLE).Details about the field campaign are described by Su et al. (2009).
The instrumentation at the forest site was similar to that at the vineyard.Net radiation was measured at a height of 35 m.Soil heat flux was measured with Hukseflux heat flux Hydrol.Earth Syst.Sci., 13,[749][750][751][752][753][754][755][756][757][758]2009 www.hydrol-earth-syst-sci.net/13/749/2009/ plates at 0.5 cm depth below the litter layer.Temperature, humidity and wind speed were measured at 35 m height.Because of an issue with the CNR1 radiometer (the temperature of the instrument was not correctly measured), the outgoing longwave radiation could not be used to estimate surface temperature.Instead, contact temperatures were measured with Negative Temperature Coefficient (NTC) sensors, 9 of which were attached to needles and branches and 8 placed on the soil.The average temperature of the 9 NTC's connected to the vegetation was used as an estimate of surface temperature.Because of the dense vegetation, the contribution of soil temperature was neglected.Contact temperature measurements were only available between 15 and 21 June 2006.Meteorological measurements were carried out and data stored at 1 min intervals.Half-hourly averages were used in this study.
Sensible and latent heat flux were measured with a CSAT3 sonic anemometer (Campbell Scientific, USA) and an open path infrared gas analyser (LI-7500, LICOR Inc., Lincoln, NE, USA) installed at 47 m height.The data were processed with the software package ECpack (http://www.met.wau.nl/projects/jep/index.html)and corrections were carried out according to Van Dijk et al. (2004).

Results
Figure 1 shows plots of a number of key variables versus time, measured and calculated with the three above mentioned methods, for one example day with clear sky conditions, for each of the three field sites.The plotted variables are the difference between surface and air temperature (T s −T a ), the aerodynamic and surface conductance, g a and g s , and the sensible and latent heat flux, H and λE.
The concuctance-approach by definition follows the a priori values for both g a and g s , without making use of the measurements of T s .For the Barrax site, this approach results in much lower modelled than measured T s −T a .The T s -approach, by definition, follows the measurements for T s , irrespective of the corresponding values for g s .The modelled vales for g s are often outside the range of values found in the literature, even negative or infinitely high.The Bayesian approach compromises between following measured values for T s and following the a priori values for the conductances g a and g s .The advantage of the Bayesian approach is that, in this way, extreme values of T s , g a or g s are avoided.
The Bayesian approach has a second advantage: the posterior values of g a and g s reveal actual information about the surface which is neither present in the a priori values, nor in the model.Posterior aerodynamic conductance g a does not deviate much from a priori, but surface conductance g s shows a clear diurnal pattern.The highest values of g s are found in the late morning for the maize and the forest.In the afternoon, g s decreases until a minimum at 19:00.During daytime hours, g s at the vineyard is lower than g s at Fig. 1. and modelled values of (T s −T a ), aerodynamic conductance g a , surface conductance g s , sensible heat flux H and latent heat flux λE, versus time.The data represent example days for a fully-grown maize field in Sonning (UK) on 13 September 2002, a vineyard in Barrax (Spain) on 15 July 2004, and a forest in Speulderbos (The Netherlands) on 16 June 2006.The thin solid line represents the T s -approach, the dashed line the conductanceapproach, the bold solid line the Bayesian approach, and the symbol "x" represents a field measurement.
the forest and the maize site.These patterns emerge even though a priori values for g s are constant and equal for all three sites.The patterns agree with conceptual understanding: semi-empirical models often predict a decreasing g s caused by stomatal closure during the afternoon.The lower g s for the vineyard can be explained by the lower vegetation cover (0.33) than that of the other two sites (closed canopies).
Figure 2 shows surface conductance g s versus vapour pressure deficit for the four sites.This figure shows a negative correlation between g s and vapour pressure deficit similar to that described in the literature (e.g.Leuning, 1995).
During the night, due to the small vapour gradient (q s −q a ), T s is relatively insensitive to g s .As a result, posterior g s returns to the a priori value during the night.This has little consequences for the calculated fluxes, since these are small during the night.
Figure 1 also shows that the conductance-approach and the Bayesian approach are unable to reproduce night time surface temperatures of the vineyard and the forest.During the night, stable, stratified air conditions are formed, in which the aerodynamic conductance strongly reduces (Massman and Lee, 2002).We do not find this strong reduction of g a in the posterior values.The reason is that the chosen standard deviation for θ 1 (σ θ1 ) is too low to include night time stable conditions.Increasing the standard deviation of θ 1 in order to account for stability is possible, but this would mean that the assumption of a Gaussian distribution of θ 1 is no longer valid.For the maize field (Sonning), for all half-hourly data in the measurement period, modelled versus measured sensible heat flux H (upper graphs) and latent heat flux λE (lower graphs), for conductance-approach (left), the T s -approach (middle) and the Bayesian approach (right).
Concerning the fluxes H and λE, the following is observed.The conductance-approach performs well during the night and in the morning, but underestimates H and overestimates λE during the afternoon.The reason is, that the a priori values for g s are constant, and that no stomatal closure is included in the conductance-approach.The T s -approach performs well during the afternoon, but poorly during the night and in the morning.The values of the fluxes calculated with the Bayesian approach always vary between that of a priori and the T s -approach, and closely follow the superior approach: the conductance-approach in the morning, and the T s -approach in the afternoon.Consequently, the Bayesian approach closely follows the measurements.Figures 3 to 5 show the results for the whole measurement period as scatter plots of modelled versus measured variables, for the three sites.For all three sites, the Bayesian approach results in the lowest root mean square error for the fluxes.The Bayesian approach reduces both the bias and the scatter compared to the other two approaches.
Both the T s -approach and the Bayesian approach overestimate λE of the vineyard.This may be caused by an eddy covariance flux measurement error rather than a model error.The vineyard was relatively small and surrounded by bare land, stubble fields and irrigated crops, which contaminated the eddy covariance signal (Timmermans et al., 2009).The high measured λE values may be attributed to an adjacent irrigated maize field.A different problem with the measurements is apparent in the forest, where measured R n was 10% higher than the sum of H , λE and G.Because the model Hydrol.Earth Syst.Sci., 13,[749][750][751][752][753][754][755][756][757][758]2009 www.hydrol-earth-syst-sci.net/13/749/2009/ forces energy balance closure (Eq.1), it is possible that λE is overestimated, while H is not underestimated (Fig. 5).
It is difficult to reproduce measured H and λE when these measurements are not used for calibration.Obviously, all of the three approaches would perform significantly better if the measured fluxes were used for calibration.However, in practice, measured fluxes are rarely available.Without measured fluxes for calibration, both the conductance-approach and the T s -approach may perform poorly (see, for example, the conductance-approach in the vineyard (Fig. 3) or the T s approach at the forest site, Fig. 5).The Bayesian always performs better than the other two approaches (a lower RMSE).Moreover, unrealistic parameter values are avoided.

Discussion and conclusion
In data assimilation, it is critical to quantify the uncertainty associated to parameters and measurements.This is not always possible for lack of data.A starting point is to identify sources of error.In this particular study, the uncertainties of aerodynamic and surface conductance are relevant.
Dimensionless "aerodynamic" parameter θ 1 can vary in space and time for two different reasons: variations in surface roughness and variations in atmospheric stability.The first is the dominant cause of variations in space (at a specific time of a satellite overpass), whereas the second is the dominant cause of variations in time (at a specific field site).Similarly, surface conductance, θ 2 , can vary in space and time for two different reasons: variations in plant species, vegetation density and soil moisture content on the one hand (spatial variations) and variations in the diurnal cycle of stomatal regulation on the other hand (temporal variations).
In this study, no distinction is made between sources of error.High standard deviations are used for the two conductances in order to cover different errors, although the effect of stability on aerodynamic conductance during the night is not incorporated (see Sect. 4).One may argue that some of the errors are of random nature, whereas others are more constant in time.An approach to handle such differences in the nature of errors is to distinguish between bias and random noise.This approach has been tested for the data of this study as well, leading to results similar to those presented in this paper.
The simple Bayesian approach led to improved estimates of sensible and latent heat flux of maize, vineyard and forest compared to more classic approaches which use either measured surface temperature or a priori parameter values.Posterior estimates reveal the diurnal pattern of surface conductance during the day (stomatal regulation), without using measured fluxes as input.

Fig. 2 .
Fig. 2. Scatter plot of posterior surface conductance versus atmospheric vapour pressure deficit (vpd) for all half-hourly data of the three sites.

Table 1 .
Parameters and variables in the energy balance and bulk transfer equations (Eqs. 1 to 3).

Table 2 .
A priori values for θ 1 and θ 2 , and standard deviations for θ 1 , θ 2 and T s .