Using a thermal-based two source energy balance model with time-differencing to estimate surface energy fluxes with day-night MODIS observations

. The Dual Temperature Difference (DTD) model, introduced by Norman et al. (2000), uses a two source energy balance modelling scheme driven by remotely sensed observations of diurnal changes in land surface temperature (LST) to estimate surface energy ﬂuxes. By using a time-differential temperature measurement as input, the approach reduces model sensitivity to errors in absolute temperature retrieval. The original formulation of the DTD required an early morning LST observation (approximately 1 h after sunrise) when surface ﬂuxes are minimal, limiting application to data provided by geostationary satellites at sub-hourly temporal resolution. The DTD model has been applied primarily during the active growth phase of agricultural crops and rangeland vegetation grasses, and has not been rigorously evaluated during senescence or in forested ecosystems. In this paper we present modiﬁcations to the DTD model that enable applications using thermal observations from polar orbiting satellites, such as Terra and Aqua, with day and night overpass times over the area of interest. This allows the application of the DTD model in high latitude regions where large viewing angles preclude the use of geostationary satellites, and also exploits the higher spatial resolution provided by polar orbiting satellites. A method for estimating nocturnal surface ﬂuxes and a scheme for estimating the fraction of green vegetation are developed and evaluated. Modiﬁcation for green vegetation fraction leads to signiﬁcantly improved estimation of the heat ﬂuxes from the vegetation canopy during senescence and in forests. When the modiﬁed DTD model is run with LST measurements acquired with the Moderate Resolution Imaging Spectroradiometer (MODIS) board


Introduction
Over the past 35 yr, a wide variety of approaches have been developed to model the surface energy balance using satellite-derived observations of land surface temperature (LST) (Kalma et al., 2008), with ongoing work in a number of techniques such as the triangle method (de Tomás et al., 2012) or one source energy balance models (Boulet et al., 2012). One of the more robust modelling approaches is the two source energy balance (TSEB) thermal-based modelling scheme, which explicitly treats the energy fluxes emanating from the soil and canopy and partitions the observed LST between the two components based on the fractional area they each occupy in the LST pixel (Norman et al., 1995). The TSEB scheme has been successfully applied for estimating surface latent and sensible heat fluxes at regional to continental scales using geostationary satellite surface radiometric temperature observations within a regional modelling system called the Atmospheric Land-EXchange Inverse (ALEXI) model (Anderson et al., 2007). The ALEXI R. Guzinski et al.: Using two source energy balance model with day-night MODIS observations modelling system addresses a critical limitation of thermalbased energy balance models regarding sensitivity to errors in absolute measurements of LST, which can be on the order of several degrees when derived from satellites due to atmospheric and surface emissivity effects. ALEXI reduces this sensitivity by using a time-differential measurement -the change of LST between two observations during the morning growth phase of the atmospheric boundary layer -which can be retrieved with better accuracy. An associated spatial disaggregation technique, DisALEXI (Norman et al., 2003), uses LST data from polar orbiting satellites to improve the spatial resolution of the modelled flux images for use in a variety of operational applications (Anderson et al., 2012a).
The Dual-Temperature Difference (DTD) model, introduced by Norman et al. (2000), also addresses the issue of sensitivity of thermal-based models to errors in absolute measurements of LST. Like ALEXI, the DTD also requires two LST observations -one early in the morning and one later in the morning or in the afternoon. However, a simpler solution scheme is employed, thereby reducing the number of required inputs and model complexity in comparison with ALEXI. The original model formulation requires an early morning LST observation (approximately 1 h after local sunrise) when fluxes are usually minimal. This means that, like ALEXI, it is dependent on the high temporal resolution of geostationary satellite measurements, which is unsuitable for applications at higher latitudes, such as in northern Eurasia and northern North America, where the view zenith angle (VZA) from geostationary satellites is large, causing loss of spatial resolution and accuracy due to longer atmospheric path lengths. These same issues also affect DisALEXI, which relies on the availability of ALEXI-derived fluxes as the normalization basis for disaggregation. The DTD has been evaluated primarily in rangelands and croplands during the growing season but before the onset of senescence. Therefore, its accuracy in forested ecosystems or other phenological stages is unknown.
In this paper, modifications to the DTD model are presented that enable it to be used with LST observations from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor aboard the Terra and Aqua polar orbiting satellites, facilitating regional surface energy flux modelling over boreal regions. First, a scheme for estimating the fraction of vegetation that is green, f g , using MODIS vegetation indices is evaluated. The green fraction is an important parameter within the model, and is used to adjust estimates of canopy transpiration based on a modified Priestley-Taylor approach . Incorporating the green fraction parameterization improves DTD model accuracy in forested ecosystems and during senescence by taking into account the phenological development of the vegetation. Next, a method for modelling the nocturnal energy fluxes is developed, taking advantage of the fact that the Terra and Aqua satellites combined provide at least two nightly and two daily acquisitions every 24 h. Finally, we consider uncertainty in flux estimates related to using the MODIS LST product. This includes adjusting the model for the different VZA associated with the day and night LST observations and considering the impact of the accuracy of the MODIS LST product.
Section 2 outlines the original DTD formulation along with the modifications proposed in this paper. In Sect. 3 we describe model validation sites, both in Denmark and in the USA, and the MODIS products used as input to the model. In Sect. 4 we evaluate the impact of the proposed modifications on the accuracy of the modelled fluxes in comparison with tower-based flux measurements in a variety of ecosystems, first running the model using in situ LST measurements and then using MODIS LST retrievals. Regional maps of energy fluxes over a hydrological observatory in western Denmark are also presented. Finally, in Sect. 5 we summarize the results and present topics for further research.

The original DTD model description
The DTD model implements the TSEB land-surface modelling scheme (Norman et al., 1995) in a time-differential mode. The directional radiometric LST, T R (θ ), is partitioned between the vegetation canopy and soil temperatures, T C and T S respectively, according to Norman et al. (1995): In Eq. (1), the fraction of view of the radiometer occupied by vegetation, f (θ ), is calculated by where θ is the VZA of the thermal sensor, F is the leaf area index (LAI) and (θ ) is the clumping factor of the vegetation at view angle θ (Kustas and Norman, 1999) and has a value of less than 1 for clumped vegetation. Using these estimates of canopy and soil temperatures, together with measurements of air temperature, T A , air density, ρ, and heat capacity of air, c p , the sensible heat fluxes for the canopy and soil, H C and H S respectively, can be derived separately. The total flux H , being the sum of the two components, can be derived by rearranging Eq. (14) from Anderson et al. (1997): In Eq. (3), R A and R S are the aerodynamic and diffusive resistances to heat transport in air above the canopy and just above the soil respectively, assuming minimal flux exchange between soil and canopy elements (i.e. flux-resistance network for soil and canopy is in parallel), and are corrected for stability using the Richardson number approximation for Obukhov length, L, and Eqs. (10) and (11) proposed in Norman et al. (2000). The sensible heat from the canopy, H C , is constrained by the canopy energy budget : where the symbol s is the slope of the saturation vapour pressure versus air temperature curve, γ is the psychrometric constant, α PT is the Priestley-Taylor coefficient with an initial value of 1.26 (Priestley and Taylor, 1972), and f g is the fraction of vegetation that is green and transpiring (see Sect. 2.2). LE C is the latent heat flux from the canopy (primarily canopy transpiration), and R n is the net radiation absorbed by the canopy calculated by Eq. (8b) from Norman et al. (2000). The net radiation of the soil and canopy system, R n , is estimated as the sum of net shortwave and long-wave radiation above the canopy. Net shortwave radiation is calculated from the measured incoming shortwave radiation and surface albedo, while net long-wave radiation is estimated from measured air temperature and LST using the Stefan Boltzman equation with atmospheric emissivity calculated as in Brutsaert (1975), and the surface emissivity estimated either from field observations or by MODIS.
In the original DTD model, Eq. (3) was applied at two times during the day: the first time approximately 1 h after sunrise, at time t 0 , and the second time later in the morning or the afternoon, at time t i . Subtracting the sensible heat flux at t 0 from that at t i gives rise to the main DTD model equation : Since the first observation is taken to be at one hour past the sunrise, the soil sensible heat flux, H S,0 = H 0 − H C,0 , is minimal and can be omitted. In practice H C,0 is also very small so the last term is also omitted, thus avoiding the need to calculate any of the fluxes or resistances at time t 0 and simplifying Eq. (5) to In the original DTD model, ground heat flux, G, at time t i was calculated as a fixed fraction (0.3) of the net radiation reaching the soil, according to Eq. (9) in Norman et al. (2000). In this study, we model G according to the scheme proposed by Santanello and Friedl (2003), which, in addition to the net radiation reaching the soil, takes into account the diurnal variation in surface radiometric temperature, assumed to be the difference between the LST observations at t i and t 0 . Finally the total latent heat flux, LE, at time t i is calculated as a residual of the other fluxes, and the latent heat flux from the soil, LE S , is also computed as a residual:

Modifications to the effective Priestley-Taylor coefficient
In combination, α PT f g form an effective Priestley-Taylor coefficient that is used to modify canopy transpiration rates computed, as in Eq. (4), from the divergence of net radiation within the canopy layer. Originally the DTD model was tested predominantly in rangelands and croplands during the growing season and before the onset of senescence . Under these conditions the model works quite well with the assumption that the vegetation is fully green (f g set to a value of 1) and transpiring at the potential rate (α PT initially set to the default value of 1.26, Priestley and Taylor, 1972). For canopies that are either stressed and not transpiring at the potential rate or are not fully green, modification to the effective α PT is required to yield reasonable partitioning between LE C and LE S . This can be accomplished by modifying either α PT or f g or both parameters, as appropriate.
In the TSEB, α PT is internally modified from its initial value if the model results in negative values in the soil evaporation rate, LE S , given that condensation on the soil is unlikely to occur during the day (Norman et al., 1995). If LE S < 0 is encountered, it is assumed that the canopy is stressed and α PT is iteratively reduced until solutions with LE S > 0 are obtained. This iterative scheme works well in ecosystems where canopy values of α PT are relatively conservative under unstressed conditions (Agam et al., 2010). For stressed canopies, the soil surface is usually dry and LE S close to zero is a reasonable assumption (Kustas and Anderson, 2009). Unless additional information about phenological condition is available, f g is typically set to unity.
In this study a number of field validation sites are located in forested ecosystems. In forests, particularly in coniferous stands, observational studies find that unstressed α PT associated with canopy is significantly lower then typical value of 1.26 (Komatsu, 2005). As a result, setting f g to unity in the DTD and downadjusting α PT iteratively from an initial value of 1.26 can lead to overestimation of LE for these ecosystems.
One approach to addressing this issue in coniferous forests is to set the initial α PT based on tree height, h C , as proposed by Komatsu (2005): Equation (9) is derived empirically and takes into account a host of physiological influences on the vegetation transpiration, including the fraction of vegetation that is green and actively transpiring. Therefore adjustment to f g = 1 is not necessary under this scheme. However this requires an estimate of h C , which is difficult to measure remotely in the absence of routine lidar datasets. Another method might assume α PT has a constant value based on an ecosystem type, which is then scaled to reflect current phenological conditions by adjusting f g . This method could be applicable not only in forests but also in grasslands and croplands during senescence. The fraction of vegetation that is green is set equal to the ratio of the fraction of photosynthetically active radiation (PAR) absorbed by the green vegetation cover and fraction of PAR intercepted by the total vegetation cover, and can be estimated using vegetation indices (VIs). Specifically, Fisher et al. (2008) proposed using the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) to estimate f g as In this study, α PT for all ecosystems was given the initial value of 1.26. Green vegetation fraction was estimated at all sites using Eq. (10).

Adapting the DTD model for night-time LST observations
By exploiting the high temporal resolution of geostationary satellites, the original DTD model was able to set the first observation time, t 0 , at one hour past sunrise when fluxes are minimal. This enabled the simplification of Eq. (5) into Eq. (6), thus avoiding the calculation of any fluxes or resistances at t 0 . Polar orbiting satellites do not have the temporal resolution of geostationary satellites. Terra and Aqua satellites used in combination provide 4 diurnal observations over most of the Earth's land surface, with Terra overpasses around 10:30 and 22:30 local time (LT) and Aqua overpasses around 13:00 and 01:00 LT. During the growing season at high latitudes, 10:30 is well past sunrise and energy fluxes are already fully developed. Therefore, the t 0 time must be associated with one of the night observations. In this study, the Aqua 01:00 LT overpass time was selected because it is closest to sunrise. At night, energy fluxes are small but they are often larger than the early morning fluxes so they could potentially influence the daytime flux estimation. During the 01:00 LT, Aqua overpass time there is no shortwave radiation component of the net radiation. Therefore the Priestley-Taylor approximation (Eq. 4) is not applicable for estimation of canopy sensible heat flux. Instead it can be assumed that at night the temperature of the canopy, T C , is close to the in-canopy air temperature, T AC , estimated by an extrapolation of the diabatic temperature profile, as described below. Then, T S can be obtained by rearranging Eq. (1) and the sensible heat of both the canopy and soil, H C and H S respectively, can be calculated from the basic TSEB equations (Norman et al., 1995): Initially T C is set to the average of T R and T A and neutral atmospheric stability conditions are assumed (|L| → ∞). This allows for calculation of resistances as described in Sect. 2.1, except in this case without the Richardson number approximation of L. Once the resistances are known, H can be calculated as the sum of H C and H S .
It is not possible to use Eq. (8b) from Norman et al. (2000) to estimate R n since it assumes solar radiation to be the dominant component of the net radiation. Instead, the longwave radiation divergences in the canopy, L N,C , and the soil, L N,S , are calculated according to Kustas and Norman (1999): where L sky , L S and L C are the long-wave emissions from the sky, soil and vegetation canopy respectively. κ L is the extinction coefficient for diffuse radiation in the canopy and is set to 0.95 for sparse vegetation (F < 1) and to 0.7 for denser vegetation (Campbell and Norman, 1998). The total net radiation at night, R n,0 , is the sum of L N,C0 and L N,S0 . L sky is calculated as described in Sect. 2.1, while L S and L C are estimated using the Stefan Boltzman equation and the modelled soil and canopy surface temperatures. At night, G cannot be estimated using the model proposed by Santanello and Friedl (2003), and instead a linear function of net radiation reaching the soil is used with the slope and intercept values of 0.3 and 35 W m −2 respectively, chosen to be similar to values found in other studies (Liebethal and Foken, 2007). LE is still calculated as the residual in the surface energy balance equation for the soil and canopy. After obtaining the first estimates of H , the night-time canopy and soil surface temperatures can be recalculated. First, the equation for diabatic temperature profile in boundary layer is used (Campbell and Norman, 1998): where 0.4 is the von Karman's constant, u * is the friction velocity, z T is the air temperature measurement height, h C is vegetation height, d 0 is the zero plane displacement length, z 0H is roughness parameter for heat transfer and H is the Monin-Obukhov stability function for heat, calculated following chapter 2.5 from Brutsaert (2005). Since d 0 ≈ 0.65h C and z 0H ≈ 0.02h C , it can be assumed that This estimate of T AC using Eq. (15) does not consider roughness sublayer effects on the temperature profile, as suggested by Harman and Finnigan (2008), but is considered a reasonable approximation for purposes of this study. u * and L can also be recalculated using the total latent and sensible heat fluxes, as described in Brutsaert (2005). The above process is repeated until both T C and L converge to stable values. If convergence is not obtained, all the fluxes are set to zero. In cases where T A is less than T R , the fluxes are also set to zero since it is not physically plausible to have unstable atmospheric conditions over land surfaces at night. Most likely in such cases there are errors in the T R retrieval.
If Eq. (5) is used with all the fluxes and resistances at time t 0 calculated as described above, then the t 0 terms on the right-hand side of the equation cancel out and the model no longer utilizes the time-differential LST observations. To avoid this, at least one of the terms with fluxes calculated at time t 0 needs to be assumed to be negligible and removed from Eq. (5). Removing these terms could potentially increase the model error when significant night fluxes are present. However, it may also increase the robustness of the model when there is bias in the temperature data, the situation that DTD was designed to address . In Sect. 4.1.2 we evaluate the impact of ignoring the night fluxes on the accuracy of the daytime estimates.

Additional considerations in using MODIS LST
When the model is driven by geostationary satellite data, the VZA remains constant between the two observations at times t 0 and t 1 . This is not the case with polar satellites, as different overpasses follow different orbital tracks and the VZA between the night and day observations changes. Therefore, Eq. (5) has to be modified slightly to take this into account: where θ 0 is the VZA of the observation at time t 0 , and θ i is the VZA of the observation at time t i .
The MODIS LST V5 products (MOD11A1 for Terra and MYD11A1 for Aqua satellites) have been validated for a number of mostly homogeneous sites using both temperature-based and radiance-based methods, and in most cases RMSE has been within 1 K Wan and Li, 2008;Coll et al., 2009).  validated MODIS LST over two lake sites and found that MODIS generally underestimated the temperature during both day and night observations and had a RMSE of 0.7 K. Wan and Li (2008) compared MODIS LST against radiance-based LST in playa, grassland, lake and bare soil sites and found again that in most cases MODIS LST underestimates ground-based measurements, both during the day and night, and errors were usually below 1 K, except for cases of bare soil where they were 1-2 K larger. Finally, Coll et al. (2009) looked at a rice field and a coniferous forest and observed an underestimation of MODIS LST of around −0.3 K with a RMSE of around 0.6 K at the rice field and negligible bias and a RMSE of 0.6 K at the forest site. The main identified causes of the underestimation bias were neglection of above-average atmospheric aerosol optical depths and difficulty in filtering out all cloudaffected observations, especially if the cloud cover was not very significant or consisted of cirrus clouds (Wan and Li, 2008). This problem is compounded at night, when it is more difficult to detect clouds (Neteler, 2010). The other identified sources of error were uncertainty in surface emissivities, especially at the bare soil and heterogeneous sites (Wan and Li, 2008) and during the periods of high soil surface surface moisture, for example after rain events (Hulley et al., 2010). No correlation, however, was found between VZA or the satellite (Aqua or Terra) and the RMSE, and in general larger errors occurred during the day than at night (Wan and Li, 2008). The prevalent bias (underestimation) in both day and night LST retrievals present in the MODIS LST product would appear to make it highly suitable for application of the DTD approach.

Danish field sites: HOBE
The Danish Hydrological ObsErvatory, HOBE, was established in 2007 to provide long-term datasets and environmental monitoring facilities for diverse hydrological studies (Jensen and Illangasekare, 2011). HOBE encompasses the Skjern River catchment on the western side of the Danish Jutland peninsula and is located in the maritime climatic zone, with mild winters and cold summers, mean annual precipitation of 990 mm and mean annual temperature of 8.2 • C. The catchment has an area of 2500 km 2 and is characterized by flat terrain, with the two main land uses being irrigated agriculture (68 % of the area) and forests (16 % of the area).
The two flux tower sites used for validation in this study are located in the Gludsted plantation forest (GLU) and Table 1. List of flux towers providing the validation datasets used in this study. Information includes the ID and name of the flux tower, the vegetation type in the vicinity of the tower, its location and the period of the data series used. Voulund agricultural site (VOU) ( Table 1). The Gludsted site is in the centre of a large homogeneous plantation dominated by 20 m tall Norway spruce (Picea abies). The eddy covariance (EC) equipment is mounted at a height of 37.5 m while the meteorological sensors are mounted 30 m above the ground. The flux measurements at GLU were collected throughout the whole of 2009, 2010 and 2011. The Voulund site is much more heterogeneous (Fig. 1). The flux tower is surrounded by crop fields, with an average area of a couple of hectares, mostly sown with winter and summer varieties of barley and maize but also with potatoes and other crops. The EC system is mounted 6 m above the ground with the flux tower footprint covering several fields, mostly west of the tower, depending on wind speed and direction and stability. The MODIS pixel also covers a number of crop fields in addition to some adjacent forest groves. The fluxes at VOU were measured from 2009 through 2011. For more detailed description of the sites and the equipment used, the reader may refer to Ringgaard et al. (2011).

ID
At both sites the fluxes and meteorological measurements were aggregated to 30 min intervals. Since DTD assumes that the incoming radiation is balanced by the outgoing energy fluxes, this balance was enforced in cases where closure was not achieved by assigning any residual energy (R n − H − LE − G) to LE, following Prueger et al. (2005). Observations where LE became negative after being assigned the residual were removed from the analysis. In addition, any observations taken on days with snow cover or very little vegetation (NDVI < 0.25) were also removed.
As a baseline case, the model was also run using predominantly ground-based measurements, with local LST estimated from the upwelling long-wave radiation measured by pyrgeometers. However, due to technical problems, no such measurements were taken at GLU and so it is not possible to run the model with local LST observations at this site. There are no field-based LAI measurements taken at the HOBE sites so MODIS LAI, smoothed using TIMESAT (Jönsson and Eklundh, 2004), was used as input into DTD.

AmeriFlux field sites
In addition to the Danish sites, a number of flux tower sites from the AmeriFlux network (http://public.ornl.gov/ ameriflux/) (Baldocchi et al., 2001) were also used to provide a more robust evaluation of model performance in different ecosystem types and climatic zones, even though these sites are at latitudes that are reasonably accommodated by geostationary satellites ( Table 1). The Black Hills (BH) tower is situated in coniferous forest, Walker Branch (WB) is in a deciduous forest, Audubon Ranch (AG) is in semi-arid grassland, Bondville (BV) is in rain-fed maize/soy bean cropland and Fort Peck (FP) is in a temperate grassland. Application period depended on archive data availability at each site. For both local tests of the DTD model (driven by towerbased measurements of LST and LAI) and experiments using MODIS-measured LST and LAI, the data from all the months of 2003 to 2006 (inclusive) have been used at BH, WB, AG sites, for 2002 to 2006 at FP, while at BV we used data from May till September for 2003 to 2006. The field measured data were 30 min averages, and energy balance closure was enforced using the residual method applied to the HOBE datasets. Measurements over bare soil or snow cover were removed from the analysis based on the same criteria as for the HOBE dataset. For more details about the AmeriFlux sites used in this analysis, see Houborg et al. (2009).

MODIS products
The modified DTD algorithm is driven mainly by a number of MODIS data products, in addition to some standard meteorological forcing (i.e. wind speed, air temperature, air pressure, relative humidity and solar radiation). Maps of LST, computed using the split-window algorithm, are distributed in the MOD11A1 (Terra satellite) and MYD11A1 (Aqua satellite) products (Wan, 2006). There are currently two versions of the product being distributed, V4 (after January 2007 V4.1) and V5, covering the period from 5 March 2000 until present-day. Among the new refinements in V5 were the relaxation of the cloudmask confidence threshold to define a pixel as clear-sky and the removal of temporal averaging to instead provide an instantaneous LST observations projected on a Sinusoidal grid . The second refinement is particularly relevant in this case since it allows for comparison of the energy fluxes with tower based measurements taken simultaneously. Although the nominal resolution of the product is 1 km by 1 km, the instantaneous field of view of the MODIS thermal-infrared sensor pixel is only that size when viewed from nadir and increases up to 2 km by 4.8 km at VZA of 60 degrees (Masuoka et al., 1998). This might be significant when evaluating modelled fluxes at heterogeneous sites. Apart from the LST layer, the M*D11A1 (MYD11A1 and MOD11A1) products contain information about the time of the observation, the VZA, and quality flag, all of which are used in the algorithm, as well as emissivity values which are needed to calculate the net long-wave radiation (see Sect. 2.1).
Another important model input parameter, used in partitioning the radiometric surface temperature and modelled fluxes between canopy and soil, is the leaf area index (LAI), which comes from the MCD15A3 product (Knyazikhi et al., 1999). MCD15A3 provides LAI estimates at 4 day temporal resolution utilizing observations from both the Terra and Aqua Satellites. The normalized difference vegetation index (NDVI) (Rouse et al., 1973) and enhanced vegetation index (EVI) (Gao et al., 2000) contained in the MOD13A2 Terra 16 day 1 km data product are used to estimate fraction of vegetation that is green as described in Sect. 2.2. Finally MCD43B3, an 8 day 1 km combined Terra and Aqua data product, is used to obtain an estimation of shortwave surface albedo, which is required to calculate the net shortwave radiation. To compute the actual surface albedo, MODISestimated black sky albedo (reflectance of direct beam radiation at solar noon) and white sky albedo (reflectance of isotropic diffuse radiation) products were combined based on the ratio of downwelling direct and diffuse shortwave radiation, which under clear skies depends on the solar zenith angle and aerosol optical depth (Jin, 2003;Lucht et al., 2000). In this application of the DTD model, an assumption is made that on clear days around the solar noon, 80 % of reflected shortwave radiation comes from a direct beam and 20 % is diffuse, which is within the observed ranges (Roderick, 1999). The effectiveness of modifying f g to improve the accuracy of modelled fluxes was tested at all the field sites with the DTD driven by field based measurements (including R n ), with the exception of NDVI and EVI used to compute f g , which were obtained from MODIS products. Since the LAI at the Amer-iFlux field sites is estimated as function of NDVI (Wilson and Meyers, 2007), it was assumed to quantify just the green part of the vegetation canopy (Houborg et al., 2009). However, when calculating properties like radiation interception or wind profile in the canopy, the whole plant area index (or total LAI) should be used. Therefore the green LAI was divided by f g to obtain plant area index in cases when f g was was less than unity; this was then used in all the DTD model equations. The same procedure was performed with the MODIS-derived LAI used at the VOU site. In the tests presented in this section, the first LST measurement was taken one hour past sunrise. Figure 2 shows the effect of adjusting f g using Eq. (10) on sensible heat flux at a coniferous forest (BH) and a deciduous forest (WB), while Fig. 3 shows the effect at an agricultural site (VOU), temperate grassland (FP) and a semi-arid grassland (AG). At all the sites the fluxes are split into green growing season and senescence phases, which is assumed to begin when green LAI decreases by 20 % after reaching its peak. At some sites there is a large variability in the dates on which the transition from growing season to senescence happens in different years due to different crop types being grown, different climatological conditions and noise present in the LAI time series. These rough estimates are sufficient for the scope of this study, in which they are mostly used for presentation and analysis purposes. In all the cases when f g was set to unity (top panels), the sensible heat flux was underestimated. At the forest sites there is no observable difference between the bias in the sensible heat fluxes during growing season and senescent periods. However the difference is very apparent at the grassland sites, with the sensible heat fluxes during senescence having a strong negative bias. At the agricultural site the senescence bias is less evident. However after the VI-based f g adjustment, the fluxes align along the 1-1 line, which is especially evident at the coniferous forest site. After f g is adjusted, the bias in H at both the coniferous and deciduous forest sites (BH and WB respectively) becomes negligible. The growing-season bias at the grassland sites (FP and AG) is minimal and the bias during senescence is reduced but still significant after the f g adjustment (Table 2). This could be due to a scale mismatch between the local LAI observations and the 1 km scale f g estimates, or the simple nature of the f g estimation scheme.
The only site where the error in the modelled sensible heat flux increases after the f g adjustment is the agricultural site, VOU, where the bias becomes positive but with a larger magnitude then when f g was kept at unity (Table 2).
To further explore the uncertainty introduced by the estimation of f g from MODIS EVI and NDVI (f g,VI ), VIbased retrievals were compared with f g derived from field observations of LAI (f g, obs ) at the BV agricultural site (Fig. 4). To compute f g, obs , it was assumed that at the beginning of the growing season, the vegetation was fully green. After LAI reached its peak, total LAI was assumed conserved and the difference between the average peak LAI and observed LAI was converted to dead leaf area: f g, obs = LAI/LAI peak (Houborg et al., 2009). Although it is not likely that LAI is conserved after reaching a peak value, since during senescence leaves shrink/shrivel or fall off, this is considered a first order approximation. During the beginning of the growing season, f g,VI shows only about 70 % of vegetation being green. At the peak of the growing season, both green fractions reach unity. The timing of the onset of senescence agrees between the field and satellite observations, but f g, obs drops to lower levels compared to f g,VI . The discrepancies between f g,VI and f g, obs have a clear effect on the estimation of the sensible heat flux at BV. Figure 5b shows that the lower values of f g,VI during the beginning of the growing season cause overestimation of H , while the higher values of f g,VI during senescence cause underestimation of H . When f g, obs values are used, the bias during both the phenological stages is minimized (Fig. 5c). This suggests that the DTD may be used to estimate fluxes during senescence in grasslands and croplands as long as changes in the green vegetation fraction are reasonably accounted for. Table 2 demonstrates the impact on bias, root mean square error (RMSE) and the coefficient of variation of RMSE (defined as RMSE divided by the mean of the observed values, CV) as a result of modifying f g at all the validation sites. Modifying f g produces, in most cases, substantial improvement in the agreement between measured and modelled H . Since the total latent heat flux is calculated as a residual, similar improvement in LE is observed, disregarding errors introduced by the calculation of G. The only exception is the increase in error at the agricultural sites during the growing season, when f g,VI is less than one even though the crops are fully green. However, the combined RMSE at BV is still significantly reduced since the magnitude of H during the senescence is much larger than during the growing season and the reduction in error during senescence is also large. At VOU the error reduces marginally during senescence but not enough to compensate for the increase during the growing season, and the overall RMSE increases slightly.
To represent optimal model performance using only remote sensing inputs, all model runs presented in the following sections used f g,VI . The only exception was at the agricultural sites (BV and VOU) during the growing season, where f g was kept at unity to better reflect the actual state of the crops.

Adjusting the model for night-time LST observations
In this section, we investigate the impact of using a first LST observation time, t 0 , in the DTD that occurs at night (01:30 LT), as in the case of MODIS applications. In these tests, nocturnal fluxes were either estimated with the model described in Sect. 2.3 or were assumed negligible as is the case in the original form of the DTD using t 0 an hour after sunrise. When the nocturnal fluxes are modelled and explicitly represented in the DTD equations, the RMSE between the estimated and observed 01:30 LT sensible heat fluxes is reduced by around 30 % for most sites (by 4-8 W m −2 at Table 2. Bias (modelled − observed) and RMSE in W m −2 and CV (unitless) of instantaneous modelled sensible heat flux at noon for fluxes during the growing season, senescence and the combined period, and using f g = 1 (f g,1 ), or f g adjusted based on VI (f g,VI ), or based on observations (f g, obs ). Model runs used t 0 at one hour after sunrise and were driven primarily by ground-based observations. MODIS products were used to estimate f g,VI at all sites and LAI at VOU.   AG, FP and VOU) and by over 50 % in the case of BH (by 23 W m −2 ). There is, however, negligible effect on RMSE at BV and an increase of almost 40 % at WB (by 10 W m −2 ), compared to the case when the nocturnal fluxes are assumed negligible.
However, using the nocturnal-flux estimates in the modelling of the day-time fluxes more often increases the discrepancies with the measurements. The 01:30 LT sensible heat fluxes are predominantly negative and taking them into account in Eq. (5) reduces the estimated daytime sensible heat flux H i . Since H i is already frequently underestimated, especially during senescence, taking the nocturnal fluxes into account actually increases the negative bias (underestimate) of the daytime H (Table 3). The above is primarily a consequence of the model, and its input data, limitations and not a reflection of any underlying physical principle. However, Table 3. Bias (modelled − observed) and RMSE in W m −2 and CV (unitless) of instantaneous modelled sensible heat flux at noon with the night-time sensible heat flux modelled (H 0,mod ) or assumed to be zero (H 0 = 0). Model runs used t 0 = 01:30 LT and were driven primarily by ground-based observations. MODIS products were used to estimate f g,VI at all sites and LAI at VOU. f g,VI = 1 during green growth at BV and VOU.   Table 3 also shows that even when nocturnal fluxes are neglected (H 0 = H C,0 = 0), the resulting accuracy is comparable to that achieved when the first DTD observation is in the early morning. This is the case for both latent and sensible heat fluxes. Therefore, application of DTD using LST data collected from the night-time MODIS overpass for t 0 does not appear to cause significant errors in the computed daytime H i .
We have also tested the robustness of the model when bias is introduced into the LST observations -a situation the DTD technique was developed to address with its use of timedifferential LST inputs. The model was run at the VOU site with both night-time and daytime LST offset by either ±1 • C or ±5 • C to reflect moderate and extreme bias conditions, and with a night offset of −1 • C and day offset of −2 • C to approximate the expected bias conditions of MODIS LST (see Sect. 2.4). Several versions of Eq. (5) were used: (i) full Eq. (5); (ii) Eq. (5) with H C,0 or H S,0 omitted depending on which one has a smaller magnitude; and (iii) both H S,0 and H C,0 omitted (Eq. 6). In all cases, night-time fluxes (H C,0 and H S,0 ) were modelled as described in Sect. 2.3. Table 4 shows that when no bias is introduced, the full form equation is the most accurate since it makes no assumption about night fluxes being insignificant; the other two versions of the equation show a small increase in the model error. When positive bias is introduced, all the versions perform very similarly and revert to Eq. (6) in case of extreme positive bias since T R,0 becomes larger than T A,0 and the night-time fluxes are set to zero (see Sect. 2.3). The biggest difference is when a negative LST offset is introduced, although the RMSE remains relatively constant for the different model versions when run with moderate LST offset. The full form of Eq. (5) experiences the largest change in bias, as expected, since the modelled night fluxes are also biased and their inclusion in the model balances the reduction in error due to the use of timedifferential LST. This is especially evident in the case of extreme negative LST offset. Equation (5) with H C,0 or H S,0 omitted performs slightly better since it neglects at least one of the biased night fluxes, while the accuracy of the model remains constant when using Eq. (6) as the day and night LST offsets cancel each other out. Similar results are obtained with the LST offset approximating the conditions of MODIS LST, with Eq. (6) being the least sensitive to the introduced offset and producing the smallest error. Therefore, it may be preferable to use this form of the model with MODIS LST data.  In summary, using LST observations at t 0 occurring during the night rather than in the early morning causes minimal degradation in the accuracy of DTD results, with some sites even showing improvement under this adjustment. This is true even when the night fluxes are assumed to be negligible, and indicates that DTD can utilize night-day observations. Therefore when using MODIS LST, the night-time fluxes can be ignored and Eq. (17) simplifies to:

Results using MODIS data
Previous studies have shown that MODIS LST estimates are often negatively biased with respect to ground-based observations, for both night and day retrievals Wan and Li, 2008;Coll et al., 2009). In this case, the use of a timedifferential LST in the first term of Eq. (18) will reduce the impact of this bias on model estimates. However, this condition is not always met at the validation sites used in this study. Specifically, when comparing the LST measured using ground instruments at the HOBE and AmeriFlux sites with the MODIS LST, the night-time MODIS LST is frequently underestimated while daytime MODIS LST is frequently overestimated as shown in Fig. 6. The biases are particularly large at VOU and WB, with FP and AG showing negligible bias during the night but strong positive bias during the day. Applying Eq. (18) in those cases compounds the biases, and may increase model errors. However, comparing MODIS derived LST to a single point-based measurement of LST is problematic for many landscapes due to the mismatch in the scale of the groundbased LST measurement and the MODIS LST pixel resolution (Coll et al., 2009). This is especially evident during the day and could be the main cause of the observed bias mismatch between the night and day measurements. Since LST is a key input for the DTD model, this can lead to a large discrepancy between the measured and modelled energy fluxes. This is further complicated by an order of magnitude difference between the source area contributing to the flux tower measurements (≈ 10 2 × 10 2 m) and the area represented by the MODIS pixel (≈ 10 3 × 10 3 m). The DTD model, modified as described in Sect. 2 and using Eq. (18), was run at the 7 validation sites using MODIS LST, LAI, f g , albedo and emissivities along with field measurements of meteorological conditions (air temperature, relative humidity, wind speed, pressure and incoming solar radiation) and vegetation height. Fluxes were modelled for a particular day only if the night Aqua LST observation was of highest quality, as indicated by the quality flag, and at least one of the day Terra/Aqua LST observations was of highest quality as well. A statistical analysis comparing the instantaneous modelled fluxes at time t i with tower measurements is presented in Table 5, with scatter plots shown in Figs. 7 (forested sites), 8 (agricultural sites) and 9 (grassland sites).  The simple R n estimation scheme in the DTD performs well at most sites, with the exception of semi-arid grassland (AG) where a large positive bias is present (Table 5). This is mostly caused by the underestimation in MODIS emissivity at the predominantly bare soil site (results not shown) and may be improved in the upcoming V6 of the M*D11A1 products (Wan, 2006). Further improvements can be expected by implementing the two-stream approach for estimating net radiation for soil and canopy . Errors in G were also acceptable (Table 5), especially within the context of the overall energy budget.
At most sites, the DTD performed well in partitioning the remaining available energy (R n − G) between sensible and latent heat fluxes (Table 5) -particularly at the forested sites (Fig. 7). Partitioning was less accurate at the grassland sites, AG and FP (Fig. 9), where strong positive bias in LE was observed, which may partly be caused by the relatively short period of vegetation activity, leading to low measured values of LE throughout most of the observation period. In addition, the overestimation of R n , particularly at AG, resulted in overestimation of LE due to its calculation as residual. The presence of irrigated crops and a river within the FP MODIS pixel likely contributed to the overestimation of LE relative to the flux tower measurements.
Two different crops were grown at BV during the four years:  Fig. 9. Instantaneous modelled net radiation (points), latent heat (diamonds), sensible heat (crosses) and ground heat (x-es) fluxes at grassland sites: AG semi-arid (a) and FP temperate (b). DTD was run using predominantly MODIS inputs with f g estimated using VI and night time fluxes ignored. 31 Fig. 8. Instantaneous modelled net radiation (points), latent heat (diamonds), sensible heat (crosses) and ground heat (x-es) fluxes at cropland sites: mixed crops (VOU, a), maize/soybean rotation (BV, b), maize/soybean with f g taken from ground observations of LAI (BV, c). DTD was run using predominantly MODIS inputs, with night-time fluxes ignored and f g in panels (a) and (b) estimated using VI during senescence and kept at unity during the growing season.   Fig. 9. Instantaneous modelled net radiation (points), latent heat (diamonds), sensible heat (crosses) and ground heat (x-es) fluxes at grassland sites: AG semi-arid (a) and FP temperate (b). DTD was run using predominantly MODIS inputs with fg estimated using VI and night time fluxes ignored. 31 Fig. 9. Instantaneous modelled net radiation (points), latent heat (diamonds), sensible heat (crosses) and ground heat (x-es) fluxes at grassland sites: semi-arid (AG, a) and temperate (FP, b). DTD was run using predominantly MODIS inputs, with f g estimated using VI and night-time fluxes ignored. the errors are reduced (Fig. 8c). In particular, the RMSE of H reduces from 101 W m −2 to 76 W m −2 and the RMSE of LE reduces from 100 W m −2 to 93 W m −2 . Nevertheless, latent heat flux is modelled quite accurately during both crop periods even when using f g,VI during senescence and f g = 1 during growing season. At the VOU site, it should be noted that the MODIS-estimated LST exhibited strong negative bias at night and positive bias during the day when compared to local ground-based measurements, which may contribute to overestimation of H and underestimation of LE by residual. The issue of inaccurate parametrization of f g is most probably present at this site as well, amplified by the different growing rates of the different ecosystems present within the MODIS pixel. Finally, the issue of sub-pixel heterogeneity is most problematic at this site -the tower is unlikely to be sampling fluxes that are representative at the 1 km scale (Fig. 1).
We further investigated the impact of using the day and night LST observations with different VZA in Eq. (18) on the accuracy of the modelled fluxes. The sensible heat fluxes modelled with MODIS input data were grouped into four categories based on θ dif = |θ 1 − θ 0 |: (1) θ dif ≤ 10, (2) 10 < θ dif ≤ 20, (3) 20 < θ dif ≤ 30, and (4) θ dif > 30. The distribution of absolute errors in each category was compared and the results are summarized in Fig. 10. The figure shows that there is no clear trend in the accuracy of the model as θ dif increases. At most of the sites there is no statistically significant difference at 5 % significance level between the medians of the absolute errors of the four categories, as indicated by the overlapping of the notched intervals. At the agricultural site BV and temperate grassland site FP, where there is significant difference between the median of the first category (with smallest θ dif ) and any of the other three categories, the first category has the smallest absolute errors. This indicates   (FP, g). The fluxes were grouped into four categories based on the absolute difference between the VZA of the day and night LST observations, θ dif : (1) θ dif ≤ 10, (2) 10 < θ dif ≤ 20, (3) 20 < θ dif ≤ 30, and (4) θ dif > 30. The boxes extend from the 25th to 75th percentiles with the central mark indicating the median. The whiskers extend to the farthest data point within 1.5 times the box range of the upper or lower box edges. Points beyond that are considered outliers and indicated as crosses. The notches on the boxes indicate the comparison interval: when the notches overlap there is no statistically significant difference at 5 % significance level between the medians. In panels (a) and (d) the comparison interval extends beyond the box due to small number of data points. that while the use of LST observations with similar VZA is preferred, this does invalidate the model output obtained with large θ dif .

HOBE area flux maps
To illustrate the utility of the modified DTD in modelling surface energy fluxes at regional scales, latent and sensible heat flux maps were produced over the Skjern River catchment within the HOBE study area (Fig. 11). MODIS inputs into the model were as described in Sect. 4.2, using LST observations acquired during a day and night Aqua overpass on the 20 April 2009. Meteorological inputs were interpolated between local measurements (Stisen et al., 2011) and vegetation height was estimated from a land cover map, with all areas classified as forest having a constant height of 20 m and other areas having height between 6 cm and 60 cm, scaled linearly using LAI. Downwelling shortwave radiation was assumed to be constant throughout the catchment and set to values measured at VOU. The modelled fluxes are instantaneous at the satellite overpass time around noon.
Although individual crop fields cannot be distinguished at the 1 km MODIS LST pixel size, and the images are early in the growing season, the maps show clear spatial patterns consistent with previous field observations ( Ringgaard et al., 2011). During 2009, the HOBE area had experienced a dry spring and summer with no significant precipitation in the 10 days preceding the 20 April. This is reflected in forests experiencing lower evapotranspiration rates compared to the irrigated crops, especially in the north and south-west of the catchment, even though they have higher net radiation due to low albedo. Pixels displaying low values of net radiation are generally located close to cloud covered pixels, and therefore likely caused by undetected clouds.
In 2011 there were 16 days during which it was possible to estimate the daytime fluxes using the DTD model over the VOU site and 23 days for GLU site. This alone is not frequent enough coverage for monitoring seasonal crop/vegetation condition or water use. However, techniques have been developed to temporally interpolate between monthly or bimonthly snapshots of ET conditions, generating full daily ET time series (Anderson et al., 2012b). Alternatively, periodic maps of the spatial patterns in water use from this method could be assimilated into soil-vegetation-atmosphere transfer (SVAT) and water balance models, providing time continuous coverage for operational applications (Crow et al., 2008). Consequently, with a few dozen satellite images and a robust interpolation technique, it may be feasible to determine seasonal water use and vegetation conditions.

Conclusions and further research
We have demonstrated the utility of the DTD model modified to take into account the issues arising from using MODIS LST data. The replacement of early morning air temperature and LST measurements with the night-time measurements as model inputs has been shown to have negligible effect on the accuracy of the modelled daytime fluxes, regardless of whether the night-time fluxes are estimated or ignored. We have therefore concluded that the nocturnal fluxes should be ignored to reduce model sensitivity to the expected MODIS LST bias conditions. We have also shown that there is no significant impact on the accuracy of the model from using LST observations with different VZAs. In addition, we incorporated a scheme for estimating the effective Priestly-Taylor α PT value by taking into account the variation in the fraction of the vegetation that is green and actively transpiring. This modification was shown to significantly improve the accuracy of the modelled fluxes at forested sites and during senescent periods at the grassland and cropland sites, even when employing a simple, VI-based parametrization. The DTD model output was evaluated at 2 flux tower sites in the HOBE observatory in Denmark and 5 AmeriFlux sites, achieving satisfactory agreement between the modelled and measured fluxes in most cases.
To resolve the issue with LST bias mismatch in the VOU site, more work should be done on validating M*D11A1 products in heterogeneous landscapes, especially croplands.
A possible alternative to validating pixel fluxes against tower measurements would be to compare catchment scale fluxes modelled with DTD versus fluxes produced by a prognostic water balance model, using independent inputs but at the same spatial scale. Additionally, a method for obtaining DTD model inputs during cloudy conditions, based for example on work by Bisht and Bras (2010), would allow for operational use of the model in environments having frequent cloud cover, such as Denmark. Finally, the robustness of method used to parametrize f g should be improved by employing a more advanced technique for estimating the phenological state of the vegetation (Tan et al., 2011).