Seasonal evaluation of the land surface scheme HTESSEL against remote sensing derived energy fluxes of the Transdanubian region in Hungary

The skill of the land surface model HTESSEL is assessed to reproduce evaporation in response to land surface characteristics and atmospheric forcing, both being spatially variable. Evaporation estimates for the 2005 growing season are inferred from satellite observations of the Western part of Hungary and compared to model outcomes. Atmospheric forcings are obtained from a hindcast run with the Regional Climate Model RACMO2. Although HTESSEL slightly underpredicts the seasonal evaporative fraction as compared to satellite estimates, the mean, 10th and 90th percentile of this variable are of the same magnitude as the satellite observations. The initial water as stored in the soil and snow layer does not have a significant effect on the statistical properties of the evaporative fraction. However, the spatial distribution of the initial soil and snow water significantly affects the spatial distribution of the calculated evaporative fraction and the models ability to reproduce evaporation correctly in low precipitation areas in the considered region. HTESSEL performs weaker in dryer areas. In Western Hungary these areas are situated in the Danube valley, which is partly covered by irrigated cropland and which also may be affected by shallow groundwater. Incorporating (lateral) groundwater flow and irrigation, processes that are not included now, may improve HTESSELs ability to predict evaporation correctly. Evaluation of the model skills using other test areas and larger evaluation periods is needed to confirm the results. Correspondence to: E. L. Wipfler (louise.wipfler@wur.nl) Based on earlier sensitivity analysis, the effect of a number of modifications to HTESSEL has been assessed. A more physically based reduction function for dry soils has been introduced, the soil depth is made variable and the effect of swallow groundwater included. However, the combined modification does not lead to a significantly improved performance of HTESSEL.


Introduction
A problem often reported in climate simulations is a systematic summer drying that results in too dry and too warm projections of summertime climate in southeastern Europe (Hagemann et al., 2004).This summer drying is associated with a strong reduction of the hydrological cycle, dry soils, strong soil evaporation and plant transpiration stress and reduced precipitation.These models often overemphasize the positive feedback between precipitation and the vapor flux due to soil evaporation and plant transpiration1 (e.g.Betts et al., 1996;Lenderink et al., 2003;Hagemann et al., 2004).Presumably, land surface processes play an important role in this feedback (Fischer et al., 2007).Improving the representation of the soil hydrological processes may impact the precipitation-evaporation feedback.Using the land surface scheme TESSEL (Tiled ECMWF Scheme for E. L. Wipfler et al.: Seasonal evaluation of the land surface scheme HTESSEL Surface Exchange over Land; Van den Hurk et al., 2000), Lenderink et al. (2003) pragmatically solved the tendency of a summer continental dry bias in their Regional Climate Model (RCM) by increasing the soil reservoir depth of TESSEL and applying a non-linear dependency of canopy resistance on available soil water.It is unclear how realistic these modifications are, and whether their application is still valid when extrapolating to changing climate conditions.
A new version of TESSEL has been developed (Hydrology-TESSEL; Balsamo et al., 2009), that appears to improve the skill of the ECMWF Integrated Forecast System to forecast the 2003 European heatwave (Weisheimer et al., 2011).A thorough test with station data and area integrated atmospheric moisture budgets (i.e.runoff data and atmospheric water balance data by Hirschi et al., 2006, andSeneviratne et al., 2004) confirmed the general improvement of HTESSEL over its predecessor (Balsamo et al., 2009).To assess the behavior of a land surface model like HTESSEL in the hydrological feedback cycle, however, a systematic evaluation of land surface evaporation at a regional scale is deemed necessary.Such an evaluation has been severely hampered by a lack of reliable and spatially explicit surface evaporation data.
To evaluate the spatial variability of seasonal mean surface evaporation from HTESSEL in a central European continental area in the Danube basin, the present study applies a spatial evaporation estimate for the 2005 growing season being derived from satellite observations.This area appears to be particularly prone to pronounced summer drying (Lenderink et al., 2003).The evaporation estimate and additional energy fluxes are derived from satellite images by using the widelyapplied, tested and validated SEBAL algorithm (Surface Energy Balance Algorithms for Land Maps) as developed by Bastiaanssen et al. (1998).
The primary objective is to assess the model skill in reproducing the spatial distribution of surface evaporation in response to a spatial distribution in precipitation forcing and a spatially heterogeneous land surface characteristics.The secondary goal is to assess the effect of a number of HTESSEL model updates.Planned simulations with the land surface scheme implemented in a full 3-D regional atmospheric model will address the model's ability to reproduce the land-atmosphere feedback in this area, which is necessary for long-term weather and climate projection.

The land surface scheme HTESSEL
In the land surface scheme HTESSEL (ECMWF, 2007;Balsamo et al., 2009) for each grid cell of the atmospheric model the land surface is represented by 6 tiles over land (bare ground, low and high vegetation, intercepted water, shaded and exposed snow).For each tile separately the energy balance is calculated: where R s and R l (W m −2 ) are the flux densities of short wave and long wave radiation, respectively, with the arrows refer to incoming (↓) and outgoing (↑) flux densities, α i is albedo, H i , λE i and G i (W m −2 ) denote the sensible, latent and soil heat flux density of tile i, respectively, λ (J kg −1 ) the specific latent heat of vaporization and E (kg m −2 s −1 ) the mass flux density of evaporation.Total H , G and λE are calculated as the area weighted average over the tiles.Soil heat is redistributed over a fixed vertical grid of 4 soil layers (extending to 2.89 m depth) using a standard diffusion scheme, allowing for thermal contributions from soil water freezing and melting (Viterbo et al., 1999).Turbulent heat and water vapor fluxes from each tile are calculated using a resistance analogy, where an aerodynamic and surface resistance accounts for the transfer efficiency of heat or water vapor over a vertical temperature and humidity gradient.The surface resistance r c is a function of R ↓ s , leaf area index LAI (m 2 m −2 ), average unfrozen soil water content θ (m 3 m −3 ), atmospheric water deficit D a (Pa), and minimum stomatal resistance r s,min (s m −1 ) (Jarvis, 1976): In particular, the sensitivity of evaporation to soil water content is relevant to discuss here, as it affects the seasonal evolution of evaporation and soil water content, i.e.: where θ wp and θ fc are the soil water contents at permanent wilting point and at field capacity, respectively, and θ is the root density weighted average water content over all soil layers of the unfrozen soil water.Hence, when θ < θ fc the resistance increases and becomes infinite at wilting point.Vertical root density distributions have been derived following Zeng et al. (1998) and adapted to a multilayer configuration.Coefficients for f 1 and f 3 are taken from a lookup table, for which an externally prescribed vegetation type forms the entry.Vegetation data are derived from ECOCLIMAP (Masson et al., 2003).
The water balance (mm d −1 ) at the land surface is described by: where W represents the change in water storage of the soil moisture and interception reservoir, S the change in accumulated snowpack, P represents precipitation, E represents evaporation of soil (E soil ), vegetation (E veg ) and intercepted water (E i ), R surface and subsurface runoff (see Fig. 1).Initially, precipitation is collected in the interception reservoir until it is saturated.Then, excess precipitation is partitioned between surface runoff and infiltration into the soil column.When the imposed water flux exceeds the maximum possible soil infiltration rate, excess water is taken as surface runoff as described by the so-called Arno scheme, while accounting for sub-grid variability related to orography (Dümenil and Todini, 1992;Van den Hurk et al., 2002).Soil water flow in HTESSEL is described by the diffusivity form of the Richards' equation using the same four-layer discretisation as for soil temperature (with increasing thickness from the soil surface downwards, i.e. 0.07 m, 0.21 m, 0.72 m and 1.89 m).The dependencies of the soil hydraulic conductivity k (m s −1 ) and soil water diffusivity D (m 2 s −1 ) on θ are described by means of the analytical functions of Van Genuchten (1980).Hydraulic coefficients are specific for six soil textures, i.e. coarse, medium, medium-fine, fine, very fine and organic.
HTESSEL does not account for either lateral exchange of soil water between the grid elements and/or irrigation.Excess water leaves the domain as either surface or subsurface runoff.At the bottom of the soil column, free drainage is assumed.Alternative lower boundary conditions are not considered.

Transdanubian test region
The test region covers the Western region of Hungary between approx.45.5-48.5 • N and 16.0-20.0• E being the Transdanubian region.Most of the area is flat and bounded by the Alps in the Southwest and the Tatra in the Northwest.The climate of Hungary can be described as a typical European continental climate with warm, dry summers and fairly cold winters.Average precipitation, P is 612 mm yr −1 and the average annual temperature at 2 m height, T a2m is about 10  19.6 • C and the average winter T a2m is 0.4  1995).About 2/3 of the land is under cultivation.The remaining vegetation is mainly deciduous forest and mixed forest (Masson et al., 2003).In Fig. 2 the percentage of areas under irrigation is given for the area considered (Siebert et al., 2007), which reveals that along the Danube valley the percentage of irrigated land is up to 50%.The annual amount of irrigated water associated with these figures is unknown, e.g. it depends on the type of crop, irrigation technique, climate and season.Measurements taken at the two flux-towers from the Car-boEuropeIP database (Tuba et al., 2005) being located in Matra and Bugac, were used as ground truth of the satellite observations.The Matra tower is located at 47 • 50 30 N and 19 • 43 33 E at 350 m a.s.l. and the tower in Bugac is located at 46 • 41 30 N and 19 • 36 06 E at 111 m a.s.l.(see Fig. 2).Both towers are situated in a grassland ecosystem.

Areal precipitation using TRMM
Area covering space-born precipitation, P , is provided by the Tropical Rainfall Measuring Mission (TRMM) on a monthly basis at a resolution of 0.25 • .Comparison of the space-born annual precipitation with the precipitation measured at 35 weather stations in the region shows that TRMM appears to overestimate the station gauge-based measurements in 2005.This occurs especially in low precipitation areas, where differences could be up to 400 mm.Therefore, we corrected the TRMM precipitation using a linear regression relation between satellite and ground observations on an annual basis (see Fig. 3a).As such we combine the spatial structure in the TRMM data with the reliable precipitation volumes of the gauge data.over the test area is given in Fig. 3b.The standard error of the corrected annual TRMM precipitation with respect to the meteostation data is about 85 mm.Annual precipitation over 2005 was on average 652 mm.The highest annual precipitation was measured in the mountainous southwestern and northeastern part of the region (around 900 mm).In the Danube valley including Lake Balaton, the annual precipitation was low, down to 450 mm.

SEBAL algorithm
The Surface Energy Balance Algorithms for Land Maps (SEBAL, Bastiaanssen et al., 1998) is applied to obtain highresolution maps of the actual evapotranspiration.SEBAL is an energy partitioning algorithm that solves the surface energy balance pixel wise from both satellite images and standard meteorological measurements.Standard inputs are firstly satellite derived maps of the Normalised Difference Vegetation Index (NDVI), surface albedo and surface temperature and secondly, measurements of the air temperature T a2m ( • C), wind speed u (m s −1 ), relative humidity RH(−) and incoming short wave radiation (R ↓ ) being obtained from meteorological stations.Updates of the algorithm's first version, as outlined in Bastiaanssen et al. (1998), have been provided by Bastiaanssen (2000) and Bastiaanssen et al. (2005).SEBAL has ever since it's development been tested and validated over a wide range of land surface types under varying climatic conditions.These included natural vegetation and agricultural crops under water-stressed as well as under well-watered conditions.A summary of these validations studies, the methods of validation, and the accuracies that were found, have been provided in Bastiaanssen et al. (2005).They found that on a daily basis the accuracy, i.e. the degree of closeness of the measurement, of the ET act estimation is around 85% whereas for seasonal or annual ET act estimations an accuracy of approximately 95% may be reached (for an accuracy of 100% the measurement and the true value are equal).A full explanation of SEBAL goes beyond the scope of this paper and reference is made to the previously cited work by Bastiaanssen.Some basic background is provided in the section below.
The latent heat fluxes are computed following Eq.( 1).Daily net radiation R n24 (W m −2 ) is computed in SEBAL using satellite measured broadband surface albedo α(−), extraterrestrial solar radiation R The a spatial resolution of 1 km.Meteorological data have been taken from 35 stations in Western Hungary and bordering countries and spatially interpolated to 1 km grids by using an interpolation method that includes land use, vegetation density and elevation (Voogt, 2006).Since R ↓ s was not available from the 35 meteorological stations, R ↓ s from the flux-towers at Bugac and Matra is used instead.The measurements from the two towers have been averaged and used as input to the SEBAL calculations, hence ignoring spatial patterns of incoming radiation.
The weekly energy fluxes were obtained by re-applying the SEBAL algorithm on a weekly basis.The weekly average meteorological data were used as input to the model, whereas the bio-physical parameters such as surface albedo, NDVI, emissivity, surface roughness and bulk surface resistance, were estimated at the time of cloud free satellite observations and assumed constant over the week.This method has been tested and validated over a wheat dominated area in Mexico (see Zwart and Bastiaanssen, 2007).

Evaluation of SEBAL daily net radiation and evaporation using tower data
For 19 cloud free observation days, daily averaged R n and λE from SEBAL have been compared with the R n and λE data obtained from the towers Matra and Bugac.At both towers, λE was measured using eddy-correlation.At both Matra and Bugac, the daily energy balance did not close, i.e. the available energy (R n +G) was larger than (H +λE).The difference on a sunny day in June could be up to 50 W m −2 .This difference is a well-known flaw of the eddy-correlation method (Wilson et al., 2002;Foken et al., 2006).The energy balance has been corrected by increasing H and λE while keeping the Bowen ratio constant.In the Fig. 4a and b the correlation is given between SEBAL and ground based corrected R n and λE, respectively.For Matra SEBAL overestimates R n by 4.8% and for Bugac by 8.4%, which is relatively good considering the difference in spatial scales between the observation methods.The difference between SEBAL and ground based λE is larger: 3% for Matra and 23.6% for Bugac, which complies with the accuracy of daily SEBAL estimates found in other studies (Bastiaanssen et al., 2005).Comparison of evaporative fraction, i.e. λE/R n , resulted in an overestimation 24% and 41% for Matra and Bugac, respectively and a coefficient of determination, R 2 of 0.14 and 0.56, respectively, which is poor for the Matra tower.

SEBAL evaporation in the Transdanubian region
In Fig. 5 the SEBAL seasonally averaged evaporative fraction λE/R n is given.Data have been downscaled to the model spatial resolution of 0.25 • .λE/R n has been derived using the total λE and R n over the growing season, which covers 30 weeks and starts in week 13 (26 March 2005).The spatial pattern of λE/R n is similar to the spatial pattern of precipitation shown in Fig. 3b, which suggests that λE/R n is to a large extent controlled by precipitation. Figure 6 shows the relationship between annual P (mm, corrected; TRMM product) and seasonal E(mm) for each grid cell.We added the calculated standard error of TRMM P (85 mm) for reference.For low P (450 < P < 750 mm), E monotonically increases with P , suggesting moisture-limited evaporation.For larger P (>750 mm) evaporation ceases to increase, pointing to radiation-controlled evaporation.The grid cells situated above the line: E = P may need additional recharge to sustain the evaporation rates (P − E) < 0. To test the statistical significance of the null-hypothesis; that P − E < 0 for the dots above the line E = P , we calculated the p-value of this null-hypothesis, using an estimated standard deviation of 85 mm.This p-value ranged between 0.52 and 0.87, which is fairly significant.
This phenomenon can also be observed in Fig. 7, which shows a map of the water balance deficit (potential recharge), i.e. the difference between TRMM P and SEBAL E. The red grid cells (P − E < 0) are situated along the river Danube, which is known to contain irrigated cropland (see Fig. 2), and could be influenced by shallow groundwater that facilitates capillary rise of water inside the soil column.These areas coincide with low precipitation areas.The blue areas, where P − E > 200 mm, are mainly characterized by mountainous terrain related to lateral (sub-)surface flow as well as lower soil thickness (and therefore reduced water availability and lower E).

Atmospheric forcing
The test domain has been divided into 170 grid cells at a resolution of 0.25 • .For this domain, simulations covering the entire year 2005 have been carried out, in which HTESSEL was forced by 3-hourly fields of precipitation P , radiation R ↓ s and R ↓ l , temperature T a2m , humidity q and wind speed at 10 m height u.These fields are taken from a simulation with the Regional Atmospheric Climate Model (RACMO2.1;Van Meijgaard et al., 2008) driven by ECMWF operational analysis.This set-up was preferred above interpolation of ERA-40 data, in order to avoid imbalances in the atmospheric driving fields originating from the data assimilation applied in ERA-40 (Uppala et al., 2005).With this set-up the right synoptic variability has been retained as well as atmospheric forcing variables that were in mutual agreement.The operational land surface scheme used in RACMO2.1 was TESSEL (Van den Hurk et al., 2000).The projected average T a2m over 2005 was 9.5 • C and the average summer T a2m was 19.8 • C.
RACMO shows systematic differences in comparison to the bias-corrected TRMM precipitation and the radiative fluxes used to drive SEBAL.To allow a meaningful comparison between HTESSEL and SEBAL mean precipitation and radiation from RACMO have been rescaled prior to running HTESSEL following a similar approach as used by Sellers et al. (1996).
The weekly averaged R n obtained through SEBAL as well as the in situ observations at Matra and Bugac are found significantly higher than the values calculated by RACMO.In Fig. 8 the weekly averaged daily R n is shown for the entire test domain.We adjusted R ↓ s obtained from RACMO for each model cell and 3-hourly timestep, such that the weekly net radiation of SEBAL and RACMO are equal: where RACMO downward and HTESSEL upward obtained with a preliminary run.Since, R ↑ l results from solving the energy balance, it is sensitive to R ↓ s .On average ξ w is larger than 1, implying the adjustment of R ↓ s to be an increase.This results in an increased R ↑ l and, hence, in a weekly R n that is generally lower than the intended R n, SEBAL w .The maximum difference is 10%.
To improve the precipitation model input we used the TRMM product (see Sect. 3.2).The monthly bias-corrected TRMM precipitation has been disaggregated using 3-hourly RACMO data in order to obtain a precipitation data series with sub-daily variation, while assuring that the areal mean is in accordance with amounts derived from the gauge-based network.The 3-hourly precipitation, P corr is: where χ m = P m,TRMM P m with P m being the monthly averaged RACMO precipitation (i.e.rainfall and snowfall) and P m,TRMM the monthly bias-corrected TRMM precipitation at the nearest data point.On a yearly basis, the scaling factor χ m ranged between 0.7 and 1.2.

Soil and vegetation data input of HTESSEL
Soil hydrologic parameters are taken from the FAO soil map and database at a spatial resolution of 5 (FAO, 1995).Soil textural information of the FAO soil types has been translated to six texture classes: coarse, medium, medium-fine, fine, very fine and organic.For each of the soil texture classes the hydraulic conductivity k and the Van Genuchten coefficients α, n and m are specified (see also Van den Hurk en Viterbo, 2003).For each grid cell the dominant soil type is used.Vegetation parameters are taken from the ECOCLIMAP vegetation map (Masson et al., 2003) at a resolution of 5 and translated to high and low vegetation tiles.

Initial conditions of HTESSEL soil state variables
Initial water in the soil system serves as a water reservoir that is available for evaporation in extended periods of low precipitation.Proper estimation of initial soil water and snow contents may have a significant effect on the modeled water balance in the succeeding season.To show this HTESSEL has been run for two sets of initial conditions.These sets contain soil moisture, intercepted water, snow water mass, snow temperature, snow density and soil temperature.Set 1 consists of initial conditions from the hindcast run of RACMO driven by ECMWF operational analysis.Set 2 uses an equilibrium initial state, obtained by cycling the model through the 2005 forcing until equilibrium was reached, i.e. using the convergence criterion of less that 1.25% difference in total soil water volume.In Table 1 mean, maximum and minimum of the initial total soil water storage (mm) and the water equivalent snow thickness (mm) in the grid cells are given for Set 1 and 2. Soil water storage differs greatly between the two sets, implying a large difference in the total annual amount of water that is available for evaporation.Set 2 has considerable higher spatial variability of initial soil water and a thicker overall snow pack than Set 1.The relatively large snow layer for Set 2 is caused by the heavy snowfall at the end of 2005.

Comparison SEBAL and HTESSEL calculations
The HTESSEL model skills to reproduce surface evaporation have been evaluated by comparing the HTESSEL evaporative fraction λE/R n with SEBAL derived λE/R n for initial www.hydrol-earth-syst-sci.net/15/1257/2011/ Hydrol.Earth Syst.Sci., 15, 1257-1271, 2011 condition Sets 1 and 2. The evaporative fractions are averaged over the growing season, starting on the 26 March 2005 and ending 30 weeks later in the same year.In Fig. 9 a flow diagram is given that summarizes the in-and output of SEBAL and HTESSEL for the evaluation.
In Table 2 the mean, variance and the 10th and 90th percentile of the seasonally averaged λE/R n are given for SEBAL and HTESSEL Sets 1 and 2. The HTESSEL mean and 90th percentile values of λE/R n correspond very well to SEBAL.The 10th percentile of HTESSEL is lower than SEBAL, which indicates a small offset towards lower λE/R n .The RMSE's of the model simulations are approximately 9% of the mean SEBAL λE/R n .
In Fig. 10a and b the difference between SEBAL and HTESSEL seasonally averaged λE/R n is given as percentage of SEBAL λE/R n for Set 1 and 2, respectively, for each grid cell.The maximum prediction error, being the difference between SEBAL and HTESSEL relative to SEBAL λE/R n is 30%.The figures indicate that initial conditions may have a considerable impact on the spatial distribution of calculated λE/R n To eliminate explicit spatial information, λE/R n values of SEBAL and HTESSEL have been ranked from low to high and subsequently plotted in Fig. 11.The figure shows that model λE/R n is slightly lower than SEBAL λE/R n for both Set 1 and Set 2. This is most prominent for low λE/R n , in spite of the large difference between the two initial condition sets.
In Fig. 12 the seasonal evaporation inferred from SEBAL and derived from the two model calculations is plotted against the annual bias-corrected TRMM precipitation.Similar to the SEBAL λE/R n, , the calculated λE/R n is precipitation dominated, especially for the initial condition Set 2. The figure further reveals that HTESSELs skill to reproduce evaporation in areas with negative potential recharge appears to be poor for both initial condition sets.

Design and evaluation of modifications to HTESSEL
To provide a rational approach to parameterization changes, Metselaar et al. (2006) analyzed the sensitivity of calculated turbulent surface fluxes to 15 different soil process parameterizations for two climates, i.e.Continental and Atlantic.The detailed and flexible soil-water-atmosphere model SWAP that is generally used for agrohydrological studies (Kroes et al., 2008), has been employed for this analysis.The analysis indicates that especially the treatment of the lower boundary condition (free drainage, irrigation, capillary rise from groundwater) and rooting depth, but also the depth of the soil column, may have a significant effect on the partitioning of radiant energy over latent, sensible and soil heat fluxes.Additionally, Metselaar et al. (2006) showed that transpiration timing strongly responds to a change of the evaporation reduction function, i.e. from a function of volumetric soil moisture content to a function of soil water pressure head.Besides, Metselaar et al. (2006) indicated that a finer mesh of the soil column yields improved convergence.
Given the results of the sensitivity analysis of Metselaar et al. (2006), we incorporated and evaluated a number of modifications to HTESSEL that are discussed below.

The effect of water stress on the canopy resistance
We changed the function f 2 in Eq. ( 3) to a (more physically based) water pressure dependent expression as: where ψ (bar) is the soil matric pressure, defined as the air pressure minus the water pressure.The matric pressure of the permanent wilting point (ψ wp ) and the field capacity (ψ fc ) is −15 bar and −0.1 bar, respectively.For ψ < ψ fc , f 2 decreases from 1 at field capacity to 0 at wilting point.In Fig. 13, the functions f 2 as defined by Eqs. ( 3) and ( 9), respectively, have been depicted as a function of ψ.
Especially in the frequently occurring ψ range between −10 and −0.2 bar the difference in reduction is large. Hydrol

Soil depth classes
To replace the fixed soil column depth of 2.89 m, spatially variable soil depths were constructed based on the Digital Soil Map of the World and Derived Soil Properties CD-ROM (Version 3.5, FAO, 1995).Given the soil type at the FAO soil map and the soil name, phase and drainage class, taxotransfer rules were used to determine the soil depth classes at the spatial resolution of the FAO map (5 ).These rules have been developed by Van Dam et al. (1994) in the framework of a European Crop Growth Monitoring System.The rationale behind these rules is that the soil depth of interest is on the one hand physically limited by rocky material below the soil column and on the other hand determined by the maximum rooting depth, which might be reduced due to rocks and/or rocky material in the soil.For example, the soil depth of lithosols is only 10 cm, and the soil depth of histosols and arenosols is 60 cm.We distinguished between five soil depth classes: soil depths of 10 cm, 60 cm, 80 cm, 100 cm and >100 cm (i.e.2.89 m).A map of soil depths in the test region is given in Fig. 14a.30% of the test region has a soil depth that is shallower than 2.89 m.For this part of the region, the original model input soil depth has to be changed to more physically realistic depths.

Shallow groundwater
Currently, HTESSEL solves the diffusivity form of the Richards' equation while assuming free drainage at the lower boundary of the soil column.As no upward flow from groundwater is possible in the current HTESSEL model, the effect of shallow groundwater is simply represented by introducing extra storage for soils with shallow groundwater.
To this end the van Genuchten retention parameter α has been changed such that the effective soil moisture at field capacity increased by 10%, i.e. the new updated situation) while retaining free drainage as the imposed bottom boundary condition.The rephrased α is solved from the expression: where θ sat, θ r are the soil moisture at full saturation and residual saturation, respectively, n and α are soil specific parameters and h fc = ψ/(ρ w g).Then: where χ = (1.1)−1/(1−1/n) .In fact, the effect of this modification is a decreased relative hydraulic conductivity for a given soil moisture content.To obtain a global map of soils influenced by shallow groundwater, the method proposed by Van Dam et al. (1994) has been applied to FAO soil type data, i.e.Gleysoils, Phaeozems, Fluvisols, Histosols, Gleyic Podsols were labeled as being groundwater affected.
Figure 14b shows the affected grid cells in the test region.

Evaluation of the HTESSEL modifications
The modified HTESSEL model skills were evaluated in a similar way as the reference HTESSEL model.Atmospheric forcing and initial conditions remained the same.Four cases have been evaluated, and compared with the reference HTESSEL.These cases reflect the modified parameterization of HTESSEL.In Case 1 the water stress function f 2 was revised.In Case 2 a variable soil depth was applied.Case 3 considers the effect of shallow groundwater.Case 4 combines Cases 1, 2 and 3. Additionally, we doubled the number of soil layers from 4 to 8 for all Cases.The number of soil layers did not affect the calculated λE/R n significantly.Based on the analysis of Metselaar et al. (2006) we expected (slightly) improved performance of the LSS.However, the diffusivity form of the Richards' equation used in HTESSEL appears to be not sensitive to the vertical discretisation.This modification is not discussed as a separate case.The evaluated cases are listed in Table 3. Statistical properties of the calculated evaporative fractions λE/R n are given in Table 4 for each case.Also the correlation coefficients between SEBAL and HTESSEL evaporative fraction are given.In Fig. 15a and b the ranked λE/R n are given for SEBAL and the cases considered for initial condition Sets 1 and 2, respectively.
For Case 1, the new water stress function resulted in a mean λE/R n and variance that are higher than the reference HTESSEL for Set 1 as well as Set 2. The increase of λE/R n , is consistent with the new f 2 function that gives less root The hydraulic properties of the soil are: k = 0.26 e −6 (m s −1 ), n = 1.25, α = 0.83 m −1 , θ max = 0.43 and θ r = 0.01.f 2 (θ ) and f 2 (ψ) have been calculated according to Eqs. ( 3) and ( 9), respectively.
water uptake reduction for similar ψ and thus increased evaporation.However, the λE/R n is larger than SEBAL λE/R n and the RMSE is therefore larger than for the reference HTESSEL.The new stress formulation does result in a new equilibrium soil hydrology regime.In order to improve both the dynamic range and the mean evaporative fraction a recalibration of other components affecting the hydrological balance is required, such as the minimum stomatal resistance and the hydraulic conductivity of the soil.For Case 2 the variable soil depth resulted in an decreased mean evaporation and a decreased 10 percentile value compared to the reference HTESSEL.The decreased evaporation is due to decreased moisture storage capacity for soil depths less than 2.89 m, i.e. in shallow soils the soil water is depleted more easily.Although additional spatial soil information is added, it only results in slightly increased variability of the calculated evaporation.
In Case 3 little effect can be observed from the additional storage to account for shallow groundwater as compared to the reference HTESSEL.The statistical properties of Case 3 are similar to that of the reference runs.We may thus conclude that the chosen parameterization does not increase the available water significantly.
For Case 4, the parameterization of the new water stress function and the additional storage due to shallow groundwater (increase of λE/R n ) are expected to balance the effect of the variable soil depth (decrease of λE/R n ).However, the effect of the new f 2 function appears to dominate the effect of the reduced soil depth.In particular this can be observed for larger λE/R n (see Fig. 15a and b and also the 90th percentiles).
For all cases the correlation between the calculated λE/R n and that of SEBAL was less than for the reference case, except for Case 2 and the initial condition Set 1 of Case 1.This outcome indicates that the model skill of HTESSEL to reproduce the spatially variable evaporation has not been significantly improved by the modifications.

General discussion and conclusions
The primary objective of this study was to assess the model skill of the land surface scheme HTESSEL to reproduce spatial patterns of surface evaporation in response to patterns in precipitation and land surface characteristics, with emphasis on the mean and spatial variability during dry (summer) periods.The secondary goal was to assess the effect of a number of model modifications.
We evaluated HTESSEL based evaporation on MODISsatellite based evaporation for the Transdanubian Region in Hungary over the growing season of 2005 (30 weeks, starting the 26 March).The energy-partitioning algorithm SEBAL has been used to calculate the energy balance terms from satellite observations.The accuracy of the SEBAL latent heat flux density λE on a seasonal basis is approx.95%, at the used spatial resolution of 0.25 • .For the land surface model, off-line atmospheric forcing variables at a 3-hourly time interval were taken from a hindcast run of the regional climate model RACMO nested in ECMWF operational analysis.In order to match with the satellite-based observations, it was found necessary to rescale downward short-wave radiation R s such that the weekly net radiation SEBAL and HTESSEL were identical in each gridcell.By doing this we eliminated the effect of net radiation as a source of spatial variability in the λE/R n patterns of SEBAL and HTESSEL.
The spatially distributed space borne precipitation was used of TRMM to rescale the forcing from RACMO.The TRMM data was corrected as it appeared to overestimate precipitation as compared to the weather stations.
The evaluation shows that, within the test region and given the available atmospheric forcing, HTESSEL predicts the seasonal energy partitioning of the incoming radiation over latent and sensible heat flux densities reasonably well for the spatial and temporal resolution and the period considered.The statistical properties of the seasonal evaporative fraction λE/R n of HTESSEL and SEBAL are of the same magnitude.However, HTESSEL slightly underestimates λE/R n , especially for grid cells with low λE/R n .The RSME of HTESSEL is 9% of the mean SEBAL λE/R n The correlation coefficient of the calculated evaporative fraction to SEBAL evaporative fraction varies between 0.58 and 0.69, depending on the initial conditions used.For individual gridcells the picture is less positive as the prediction error is up to 30% of the SEBAL λE/R n .
The results above are based on an atmospheric forcing of which the accuracy is unknown.Since the effect of net radiation on the spatial pattern is eliminated, especially precipitation P may have a large impact on calculated evaporation E pattern (see also Fig. 12): a slight change in precipitation may change the calculated evaporative fractions.It is recommended to reevaluate the use of TRMM data as a source of spatially distributed precipitation.In addition, longer evaluation periods are needed to confirm these observations.Next to radiation and precipitation, a third source of spatial variability is the initial water in the terrestrial system available for evaporation, i.e soil moisture distribution and snow cover.As these conditions were unknown we used 2 initial condition sets to assess the effect of differences in spatial distribution on the model performance.
These sets differ in the initial conditions of soil water and snow cover taken from a RACMO hindcast run in Set 1 and an equilibrium state over 2005 in Set 2. The use of the two different initial condition sets, resulted in different spatial distributions of λE/R n .A closer look shows a great similarity between the spatial pattern of the potential recharge (Fig. 7) and the relative prediction error (accuracy) using Set 2 (Fig. 10b).This finding is confirmed by the calculated correlation coefficients between potential recharge and relative prediction error, which are 0.55 and 0.8 for Set 1 and 2, respectively.HTESSEL's ability to predict in areas where P − E < 0 appears to be lower than in areas where P − E > 0. During drier years than 2005 (which had an annual precipitation anomaly of 40 mm compared to the climatological mean of 612 mm), the prediction error may become larger.
Figure 12 illustrates the importance of representing initial terrestrial water storages correctly for modulating wet and dry meteorological anomalies.The figure shows yearly precipitation P plotted against seasonal evaporation E of each individual grid cell.E of Set 1 has a scattered relationship to P , whereas Set 2 shows a largely linear relationship.Due to the cycling over 2005, the initial state of Set 2 reflects only the signature of the atmospheric forcing over 2005, which is dominated by P in the region considered.Instead, in Set 1 the initial water stored in the soil and snow pack reflects the signal of longer-term meteorological conditions, which is probably more realistic.In grid cells with low precipitation over 2005, the effect on λE/R n may be moderated by the relatively wet soil moisture conditions originating from a previous (winter) period.The overall effect is a more scattered relationship between P en E.
The simulations performed with HTESSEL reveal a relatively low ability of the model to correctly predict E in areas where P − E < 0. However, longer evaluation periods are needed to confirm the results and the effect of the uncertainty in the precipitation patterns used needs to be evaluated.Since low precipitation areas coincide with irrigated areas, the underestimation of λE/R n in grid cells with (P − E) < 0 might also point towards enhanced evaporation due to irrigated cropland.HTESSEL does neither incorporate the effect of irrigation, nor the effect of shallow groundwater on the soil water balance.Especially, during periods with high temperatures and low humidity, additional evaporation is expected due to the availability of irrigation and groundwater.It is recommended to further assess the effect of irrigation on the land surface performance as taking into account of irrigation may improve the performance.
Like many other Land Surface Schemes (LSS), HTESSEL does not allow for lateral redistribution of precipitation by surface and subsurface flow.LSSs are now being modified to include lateral flow, groundwater flow and surface water (e.g.Fan et al., 2007;Miguez-Macho et al., 2007).These developments however are still in an experimental stage, which is largely due to difficulties in obtaining the required hydrological data that needs global coverage and a correct resolution.By bridging the gap between hydrological E. L. Wipfler et al.: Seasonal evaluation of the land surface scheme HTESSEL and climate models and thus by incorporating lateral flow, including groundwater flow, irrigation and river routing, LSS skills may significantly improve.
Based on earlier sensitivity analysis of soil hydrologic processes (Metselaar et al., 2006) we modified some of the model concepts and we added spatial information on soil depth and groundwater depth to improve the model performance, i.e. (i) revised the parameterization of the reduction of evaporation for dry vegetation, (ii) replaced the fixed soil depth with more realistic and variable soil depths, and (iii) introduced additional water availability due to capillary rise from shallow water tables.These modifications were evaluated again against SEBAL derived λE/R n using the same atmospheric forcing and initial conditions.The evaluation did not show to significant performance improvement for both initial condition sets: Modification (i) increased λE/R n too much, especially for grid cells in the higher λE/R n range.(ii) decreased λE/R n (especially in the lower range) and increased the RMSE.The unrealistically large soil thickness in HTESSEL seems to compensate for the strong reduction of root water uptake under dry conditions.At higher spatial resolutions the spatial variability of soil and vegetation characteristics may become more important and hence a more physically based description of soil moisture movement may be warranted.(iii) did not lead to significant changes in λE/R n .For the implementation into the RCM RACMO, a more rigorous parameterization for the groundwater dynamics will be needed to improve the models ability of predicting evaporative fractions in regions affected by shallow groundwater.
The evaluation of the performance of the modified HTESSEL as presented in this paper should be confirmed using a longer evaluation period.In addition, the uncertainty in the precipitation forcing of HTESSEL needs to be assessed.

↓
s,exo (Wm −2 ) and incoming short wave radiation (ground truth) R ↓ s according to De Bruin and Stricker (2000): Fig. 4. (a) Correlation between groundbased daily R n and R n -SEBAL.(b) Correlation between groundbased corrected E and E-SEBAL.The groundbased energy fluxes have been obtained from the meteorological towers Matra and Bugac, of which the energy balances have been closed proportional to the Bowen ratio.The SEBAL observations have a spatial resolution of 1 km.

Fig. 5 .
Fig. 5. Map of seasonally averaged λE/R n for the Transdanubian test region in 2005.λE/R n has been derived from satellite images using SEBAL.The spatial resolution is 0.25 • .The season extends over 30 weeks, starting in week 13.

Fig. 6 .
Fig. 6.Seasonal SEBAL derived evaporation E (mm) over a 30 week period in 2005 starting at week 13 and ending at week 43 and annual TRMM precipitation P (mm) over 2005 for each grid cell in the test region.Each point represents a grid cell.The line E = P is given for reference as well as the dotted lines representing E = P ± the standard error of corrected TRMM P .

Fig. 7 .
Fig. 7. Potential recharge (Annual P TRMM -Seasonal E SEBAL ) over 2005 in the Transdanubian region.The areas A and B represent irrigated cropland area and C represents lake Balaton.The spatial resolution is 0.25 • .The SEBAL evaporation is averaged over 30 weeks starting in week 13.

Fig. 8 .
Fig. 8. Weekly averaged, daily net radiation R n over 2005, as obtained by the SEBAL algorithm and from the HTESSEL/RACMO simulation, respectively, being averaged over the test area.The meteorological forcing and initial conditions were obtained from a one-year RACMO2-hindcast run driven by ECMWF operational analyses.

Fig. 9 .
Fig. 9. Flow diagram of the in-and output of SEBAL and HTESSEL.

Fig. 10 .
Fig. 10.Difference between SEBAL and HTESSEL seasonally averaged evaporative fraction, λE/R n as percentage of SEBAL λE/R n for initial condition Set 1 and 2, respectively.λE/R n has been derived by using the λE and R n over a 30 period starting at week 13 and ending at week 43, 2005.The blue cells refer to λE/R n overprediction, the red cells to underprediction by the model.

Fig. 11 .
Fig. 11.HTESSEL seasonally averaged evaporative fraction λE/R n values for each grid cell being ranked from low to high values and plotted against ranked SEBAL λE/R n .The SEBAL evaporation is averaged over 30 weeks starting in week 13.

Fig. 12 .Fig. 13 .
Fig. 12. Seasonal SEBAL and HTESSEL derived E (mm) over a 30 week period starting at week 13 and ending at week 43, 2005 and annual TRMM precipitation (P , mm) over 2005 for each grid cell in the test region.Each symbol represents one grid cell.The dotted line E = P and SEBAL E are given for reference.Like Fig. 6, the figure shows a correlation between P and E, especially for Set 2. For both initial condition sets the transition between precipitation dominated E and radiation dominated E is clearly visible.
Fig. 14.(a) Soil depth classes and (b) groundwater affected soils (dark cells) occurring in the test region, as based on FAO soil classification and expected rooting depth being aggregated to the spatial resolution of 0.25 • .
• C (Szalia et al.,  2005).The soils in the area can be classified as acid and nonacid loamy, well-drained soils, salt affected, sodium rich and imperfectly drained soils (Dobris report, Soil map of Europe,

Table 1 .
Mean, minimum and maximum model grid values of initial soil water storage (mm) and initial water equivalent snow (mm) within the test region.Set 1 is the initial condition set that originates from the RACMO hindcast run, Set 2 is the equilibrium initial condition state.

Table 2 .
Mean, 10th and 90th percentile of seasonally averaged λE/R n from SEBAL and HTESSEL for both the initial condition Set 1 and Set 2 using a spatial resolution of 25 km.The RMSE of the HTESSEL model predictions is given in the last row.

Table 3 .
Evaluated combinations of the modifications to HTESSEL.Four configurations (cases) were considered.The differences with respect to the reference HTESSEL is indicated in grey.f 2 dependency refers to the function used to calculate the effect of water stress.Soil depth refers to variable soil depth, groundwater effect refers to the considered effect of shallow groundwater and the number of compartments refers to the vertical discretisation of the soil column in HTESSEL.

Table 4 .
Summary statistics of the evaluated cases.
Seasonally averaged ranked SEBAL λE/R n plotted against ranked HTESSEL λE/R n for all grid cells in the test region for: (a) initial condition Set 1 that refers to an initial soil state variable condition set and (b) initial condition Set 2 that refers to equilibrium of state variables over 2005.λE/R n has been calculated by using the λE and R n over 30 weeks starting at week 13.