The Tiled ECMWF Scheme for Surface Exchanges over Land (TESSEL) is used operationally in the Integrated Forecast System (IFS) for describing the evolution of soil, vegetation, and snow over the continents at diverse spatial resolutions. A revised land surface hydrology (H-TESSEL) is introduced in the ECMWF operational model to address shortcomings of the land surface scheme, specifically the lack of surface runoff and the choice of a global uniform soil texture. New infiltration and runoff schemes are introduced with a dependency on the soil texture and standard deviation of orography. A set of experiments in stand-alone mode is used to assess the improved prediction of soil moisture at the local scale against field site observations. Comparison with basin-scale water balance (BSWB) and Global Runoff Data Centre (GRDC) datasets indicates a consistently larger dynamical range of land water mass over large continental areas and an improved prediction of river runoff, while the effect on atmospheric fluxes is fairly small. Finally, the ECMWF data assimilation and prediction systems are used to verify the effect on surface and near-surface quantities in the atmospheric-coupled mode. A midlatitude error reduction is seen both in soil moisture and in 2-m temperature.
A correct representation of the soil water buffering in land surface schemes used for weather and climate prediction is essential to accurately simulate surface water fluxes toward both the atmosphere and rivers (van den Hurk et al. 2005; Hirschi et al. 2006a). Moreover, the energy partition at the surface is largely driven by the soil moisture, which directly influences the Bowen ratio.
Differences in the treatment of the surface energy balance and soil hydrology between different land surface schemes clearly emerged in several offline land surface model intercomparison experiments, often carried out under the umbrella of the Project for Intercomparison of Land-Surface Parameterization Schemes (PILPS; Henderson-Sellers and Dickinson 1992; Henderson-Sellers et al. 1995). Also the European Centre for Medium-Range Weather Forecasts’ (ECMWF) Tiled ECMWF Scheme for Surface Exchanges over Land (TESSEL; van den Hurk et al. 2000) participated in a suite of these PILPS-type experiments, in particular the Thorne–Kalix experiment (Nijssen et al. 2003) and the Rhone Aggregation experiment (Rhone-AGG; Boone et al. 2004). In these experiments, a couple of weak components of TESSEL became evident and were further explored by van den Hurk and Viterbo (2003). Specifically, the single global soil texture, which does not characterize different soil moisture regimes, and the Hortonian runoff scheme, which produces hardly any surface runoff, were identified as priority development areas.
The effects of allowing variable soil textures and the benefit of alternative soil hydraulic parameterizations instead of the often used Clapp and Hornberger (1978) set of equations have been demonstrated in a number of studies (e.g., Kleidon and Heimann 1998; Shao and Irannejad 1999). Many land surface models (LSMs) of today’s climate and NWP models use a variable soil and/or rooting depth [e.g., Organizing Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE), de Rosnay et al. 2002; Interactions between Soil, Biosphere, and Atmosphere (ISBA), Noilhan and Planton 1989; Variable Infiltration Capacity (VIC), Wood et al. 1992; Noah, Chen and Dudhia 2001; and Community Land Model (CLM), Oleson et al. 2008]. In NWP applications (the primary focus of ECMWF’s core activity), it has never been clearly demonstrated that the use of a global constant soil texture puts limitations, for instance, on near-surface temperature or humidity forecasting skill. However, application of TESSEL outside the domain of short-term weather forecasting puts stronger demands on the ability to adequately describe the actual soil water storage, and the exploration of this variable soil texture is definitely worthwhile.
Similarly, the TESSEL approach toward surface runoff treatment shows obvious shortcomings in both the seasonal mean magnitude of runoff fraction (runoff divided by total precipitation plus snowmelt) and the temporal variability at short (1–20 days) time scales. Decharme and Douville (2007) presented the different merits of orographic-based subgrid runoff parameterizations, indicating the variable infiltration capacity approach is suitable for global hydrological modeling. Van den Hurk and Viterbo (2003) demonstrated the benefit of including a variable infiltration capacity approach for the Thorne–Kalix experiment.
Although a repetition of this experiment for Rhone-AGG actually showed that this modification yielded worse scores in terms of Nash–Sutcliffe discharge efficiency coefficients, the time spectrum of the simulated discharge in the new model version was in much better agreement with the observations than the traditional TESSEL approach. This approach has been implemented in many present-day land surface models already [VIC, Wood et al. 1992; ISBA, Noilhan and Planton 1989; Max Planck Institute (MPI), Dümenil and Todini 1992; Hagemann and Gates 2004] and is considered to describe the effect of subgrid variability on the runoff characteristics much better than the uniform infiltration saturation in use in TESSEL.
Thus, a revised formulation of the soil hydrological conductivity and diffusivity, spatially variable according to a global soil texture map, and surface runoff based on the variable infiltration capacity approach are the proposed remedies. This paper discusses the implementation of these components into TESSEL and performs an extensive evaluation in a range of experimental designs covering various temporal and spatial scales.
Offline (or standalone) verification is a convenient framework for isolating the benefits of a given land surface parameterization. A set of field site experiments and two land surface intercomparison experiments over large domains are considered. A Sahelian site and a boreal forest site have been chosen to show relevant effects of the new hydrology. Two major land surface intercomparison experiments—the second Global Soil Wetness Project (GSWP-2) (Dirmeyer et al. 1999, 2002; Gao et al. 2004) and the Rhone-AGG project (Boone et al. 2004)—provided spatialized near-surface forcing for land surface models and have been rerun with the new scheme to evaluate the water budget for accumulated quantities. In the GSWP-2 simulations, both terrestrial water storage estimates and the river discharge are examined for a number of basins. Hydrological consistency on the monthly time scale is verified. The Rhone-AGG simulation is used to examine the fast component of runoff at the daily time scale.
The coupling between the land surface and the atmosphere is then also evaluated. This is an essential step to confirm the results obtained offline, as indicated by the Atmospheric Model Intercomparison Project (AMIP) experience and its comparison to PILPS experiments (Qu and Henderson-Sellers 1998). Moreover, Koster et al. (2004) indicated a strong coupling between soil moisture and precipitation over large continental areas, generalizing the studies of Beljaars et al. (1996); therefore, soil hydrology modifications may affect the overall model climate.
To assess the impact of the new parameterization, a set of long-term atmospheric coupled integrations (13 months) with specified sea surface temperature is produced. This configuration, named climate simulation, allows for evaluating surface–atmosphere feedbacks and for focusing on the effect of the land surface modification. Annual and seasonal averages are compared to a number of independent datasets, with a focus on boreal summer months when a larger effect of the soil hydrology is expected. Finally, since in the NWP application the coupled system is subject to cyclic correction by data assimilation, an overall assessment is provided by the comparison of the land surface analysis increments in the old and new version of the land surface model. A reduction of increments between the two model versions can be interpreted as an overall improvement of the land surface representation. In this case, the soil moisture increments are considered in a long (seven months) data assimilation experiment.
In section 2, the hydrology of the TESSEL land surface scheme and the hydrology TESSEL (H-TESSEL) revision are illustrated. Sensitivity experiments are conducted to show the effect of the new parameterization on surface runoff and soil water transfer. Section 3 evaluates the soil moisture range associated with the new physiographic values (for the permanent wilting point and the soil field capacity) for a number of sites and presents the validation results at two contrasting field sites that illustrate the main behavior of the H-TESSEL and TESSEL schemes. In section 4, the regional-to-global offline simulations are introduced, together with the main verification datasets provided by the basin-scale water balance (BSWB; Seneviratne et al. 2004) the Global Runoff Data Centre (GRDC) (Fekete et al. 2000) and the 40-yr ECMWF Re-Analysis (ERA-40; Uppala et al. 2005) datasets.
Results of atmospheric-coupled simulations and data assimilation experiments are presented and discussed in section 5, together with the relevant lessons learned from the ERA-40 reanalysis. A summary of the changes and the conclusions are provided in the last section, while the datasets for the global implementation of the revised hydrology scheme are presented in the appendix.
2. TESSEL hydrology
The TESSEL scheme is shown schematically in Fig. 1a. Up to six tiles are present over land (bare ground, low and high vegetation, intercepted water, and shaded and exposed snow) and two over water (open and frozen water), with separate energy and water balances.
The vertical discretization considers a four-layer soil that can be covered by a single layer of snow. The depths of the soil layers are in an approximate geometric relation (7 cm for the top layer and then 21, 72, and 189 cm), as suggested in Deardorff (1978). Warrilow et al. (1986) have shown that four layers provide a reasonable compromise between computational cost and the ability to represent all time scales between one day and one year. The soil heat budget follows a Fourier diffusion law, modified to take into account soil water freezing/melting according to Viterbo et al. (1999). The energy equation is solved with a net ground heat flux as the top boundary condition and a zero flux at the bottom. An interception layer accumulates precipitation until it is saturated, and the remaining precipitation (throughfall) is partitioned between surface runoff and infiltration. Subsurface water fluxes are determined by Darcy’s law, used in a soil water equation solved with a four-layer discretization shared with the heat budget equation. The top boundary condition is infiltration plus surface evaporation, free drainage is assumed at the bottom, and each layer has an additional sink of water in the form of root extraction over vegetated areas.
In each grid box, two vegetation types are present: a high and a low vegetation type. An external climate database is used to obtain the vegetation characteristics, based on Global Land Cover Characteristics (GLCC) data (Loveland et al. 2000; available online at http://edcsns17.cr.usgs.gov/glcc/). The nominal resolution is 1 km. The data provides for each pixel a biome classification based on the Biosphere–Atmosphere Transfer Scheme (BATS) model (Dickinson et al. 1993), and four parameters have been derived for each grid box: dominant vegetation type (TH and TL) and the area fraction (AH and AL) for each of the high- and low-vegetation components, respectively.
where ρw is the water density (kg m−3), Fw is the water flux in the soil (positive downward; kg m−2 s−1), and Sθ is a volumetric sink term associated to root uptake (m3 m−3 s−1), which depends on the surface energy balance and the root profile (Viterbo and Beljaars 1995).
The liquid water flow, Fw, obeys Darcy’s law, written as
where λ (m2 s−1) and γ (m s−1) are the hydraulic diffusivity and hydraulic conductivity, respectively.
The top boundary condition is given by precipitation plus snowmelt minus bare ground evaporation minus surface runoff. The bottom boundary condition assumes free drainage. Abramopoulos et al. (1988) specified free drainage or no drainage, depending on a comparison of a specified geographical distribution of bedrock depth, with a model-derived water table depth. For the sake of simplicity, the assumption of no bedrock everywhere has been adopted.
TESSEL adopts the Clapp and Hornberger (1978) formulation of hydraulic conductivity and diffusivity as a function of soil water content [see also Mahrt and Pan (1984) for a comparison of several formulations and Cosby et al. (1984) for further analysis]:
where bc is a nondimensional exponent, and γsat and ψsat are the values of the hydraulic conductivity and matric potential at saturation, respectively. A minimum value is assumed for λ and γ corresponding to permanent wilting-point water content.
Cosby et al. (1984) tabulate best estimates of bc, γsat, ψsat, and θsat for the 11 soil classes of the U.S. Department of Agriculture (USDA) soil classification, based on measurements over large samples. Viterbo and Beljaars (1995) adopted an averaging procedure to calculate for a medium-textured (loamy) soil used in TESSEL the values of γsat = 0.57 × 10−6 m s−1, bc = 6.04, and ψsat = −0.338 m, compatible with the Clapp and Hornberger expression for the matric potential
The water transport in frozen soil is limited in the case of a partially frozen soil, by considering the effective hydraulic conductivity and diffusivity to be a weighted average of the values for total soil water and a very small value [for convenience, taken as the value of Eq. (4) at the permanent wilting point] for frozen water as detailed in Viterbo et al. (1999). The soil properties, as defined above, also imply a maximum infiltration rate at the surface, defined by the maximum downward diffusion from a saturated surface. In general, when the water flux at the surface exceeds the maximum infiltration rate, the excess water is put into surface runoff. The general formulation of surface runoff can be written as
where Imax is the maximum infiltration rate, T the throughfall precipitation, and M the snowmelt. Different runoff schemes differ in the formulation of the infiltration. The maximum infiltration or Hortonian runoff represents the runoff process at local scales. In TESSEL, the maximum infiltration rate Imax is calculated as
where ρw is the water density and z1 is the depth of the first soil model layer (7 cm). At typical NWP model resolutions, this scheme is active only in the presence of frozen soil, when downward soil water transfer is inhibited; otherwise, it hardly ever produces surface runoff, as shown in Boone et al. (2004).
a. The H-TESSEL revision
The H-TESSEL scheme includes the following revisions to the soil hydrology: (i) a spatially varying soil type replacing the single loamy soil; (ii) the van Genuchten (VG) formulation of soil hydraulic properties replacing the Clapp and Hornberger (CH) scheme; and (iii) the surface runoff generation changing according to a variable infiltration capacity based on soil type and local topography.
In Fig. 1b, the H-TESSEL changes are illustrated: in two adjacent model grid points with the same land surface conditions and receiving an equal amount of precipitation, the surface runoff will be different and proportional to the terrain complexity, while the soil water drainage will depend on the soil texture class.
The van Genuchten (1980) formulation provides a closed-form analytical expression for the conductivity, given as a function of the pressure head h, as
where α, n, and l are soil texture–dependent parameters. Pressure head h is linked to the soil moisture by the expression
The VG scheme is recognized among soil physicists as capable of reproducing both the soil water retention and the hydraulic conductivity and has shown good agreement with observations in intercomparison studies (Shao and Irannejad 1999). Table 1 lists parameter values for six soil textures for the VG scheme. HTESSEL uses the dominant soil texture class for each grid point. This information is taken from the Food and Agriculture Organization (FAO) dataset (FAO 2003), as detailed in the appendix, provided at a nominal resolution of about 10 km. The permanent wilting point and the soil field capacity are obtained by a specified matric potential of ψ(θpwp) = −15 bar and ψ(θcap) = −0.10 bar, respectively.
In Table 2, the volumetric soil moistures associated with each soil class are shown for saturation, field capacity, and wilting point. Also shown are the plant available water content and the percentage of land points in each class. The last row shows the corresponding values for the single loamy soil used in the CH formulation in TESSEL. Note that the plant available soil water is greater for all the new soil classes in H-TESSEL. Figure 2 shows the soil hydraulic diffusivity and conductivity for the TESSEL CH formulation and the six VG soil texture classes in H-TESSEL. In TESSEL, those were not allowed to fall below their wilting point values. At saturation, TESSEL has the highest diffusivity and conductivity. The reduced values for fine soils in HTESSEL reduce the infiltration of water and consequently the baseflow.
A variable infiltration rate, first introduced in the so-called Arno scheme by Dümenil and Todini (1992), accounts for the subgrid variability related to orography and considers that the runoff can (for any precipitation amount and soil condition) occur on a fraction s of the gridpoint area S,
where W and Wsat are vertically integrated soil water contents (θ and θsat) over the first 50 cm of soil, defined as an effective depth for surface runoff. Parameter b is spacially variable, depends on standard deviation of orography (σor), and is allowed to vary between 0.01 and 0.5. The parameters σmin and σmax are set to 100 and 1000 m, respectively, as in van den Hurk and Viterbo (2003).
The surface runoff is obtained by the Hortonian runoff formulation by integrating Eq. (10) over the grid box:
Whenever rain or snowmelt occurs, a fraction of the water is removed as surface runoff. The ratio runoff/precipitation scales with the standard deviation of orography and therefore depends on the complexity represented in the grid box as well as on soil texture and soil water content via W and Wsat. In Fig. 3, the response to a 10 mm h−1 rain rate for the six VG soil types and for the CH case in TESSEL is shown as a function of the b parameter. At field capacity, the surface runoff may vary from roughly 1%–50% of the rainfall (snow melting) rate, generally increasing with finer textures and orographic complexity.
3. Field site experiments
a. Soil properties
To evaluate the specified values for the permanent wilting point and field capacity used in H-TESSEL, we considered field observations of an agrometeorological network, a subset of the Global Soil Moisture Data Bank (Robock et al. 2000). These parameters are in fact crucial to capturing the seasonal and synoptic variability of soil moisture. The agrometeorologic networks operated by Russia (63 stations) and Ukraine (96 stations), as plotted in Fig. 4, are used for a field site–based verification of the physiographic properties. Measurements of physical soil properties are made at depths of 20 and 100 cm and comprise the volumetric density, total water holding capacity, field capacity, and level of wilting. Vegetation cover includes maize, winter wheat, and spring wheat fields. Measurements are taken periodically during the agricultural season, which starts at the beginning of the field work (April) and lasts until harvest time. In fields with a spatially inhomogeneous soil structure, several cross sections are made to obtain a reliable estimate of these quantities.
Figures 5 and 6 show the histogram of differences between the wilting level, the field capacity, and the water holding capacity (i.e., the difference between field capacity and wilting level) based on the field observations at 20-cm depth. The first panel (top) shows the histograms for the FAO soil texture–derived values for the field capacity calculated by setting ψ(θcap) = −0.10 bar. The second panel (middle) shows the same data, but the field capacity was set to ψ(θcap) = −0.33 bar. The third (bottom) panel shows the histograms for the uniform soil parameters used in the TESSEL scheme. The bias observed in the wilting point of the FAO soil map is clearly substantially smaller than the bias of the uniform value used for the wilting point in TESSEL. For both networks, the bias decreases by 50% or more when considering H-TESSEL, with values 3.0% by volume (vol%) and 1.1 vol% for Ukraine and Russia, respectively (for TESSEL, the bias value for each is 6.2 vol% and 3.4 vol%, respectively). For the field capacity, the results depend largely on the setting of ψ(θcap). Setting the value to −0.10 bar leads to differences similar to those observed for TESSEL. For the Ukrainian data, the absolute bias slightly increases by 1.3 vol%; for Russia, the absolute bias decreases by 1.2 vol%. Although Hillel (1982) indicates a matric potential ψ(θcap) of −0.33 bar as a common value for medium-textured soil field capacity, this leads in H-TESSEL to a substantial degradation when compared to the field capacity estimates. The bias increases to a high value of −4.7 vol% and −5.4 vol% for Ukraine and Russia, respectively. These biases are significantly larger than those observed for the uniform TESSEL value. Therefore, the value of −0.10 bar used for coarse-textured soil is adopted for all the soil texture classes. The water holding capacity reflects the positive effect of this choice and of the new soil scheme. For both networks, the bias decreases significantly for the H-TESSEL scheme.
b. Soil moisture and fluxes
The validation focuses on well-instrumented sites where various types of observations are available and where the soil type was likely to play a significant role in the surface model behavior. Long time series (one year or more) of hourly observations of low-level wind, temperature, and moisture, together with precipitation and downward radiation, are used as input to the land surface. Two different sites are considered for offline simulations using both TESSEL and H-TESSEL. Table 3 reports the location and the physiography (soil and vegetation) for each site. These quantities were assigned by the operational high-resolution forecast model (at about 25 km globally), adopting the nearest grid box. The sites correspond to different soil textures and to extremely contrasting climatic conditions. The initial conditions for the land surface state, including soil moisture, have been obtained via perpetual-year integrations until the schemes reach equilibrium.
1) SEBEX Sahel
The Sahelian Energy Balance Experiment (SEBEX; Wallace et al. 1991) is characterized by a desert climate with fallow savannah vegetation and sandy soil. Volumetric soil moisture reaches extremely low values in the dry season. This dynamical range is typical for sandy soil and cannot be represented by the medium texture used in TESSEL, resulting in a much wetter range. Observations of both soil moisture and evaporation are available for verification, as previously considered in van den Hurk et al. (2000). Figure 7 shows the time series of soil moisture. Measurements were taken weekly at several depths and here are shown at 17 cm below the surface and thus compared to the soil model layer 2 representative of the depth range 7–28 cm. It is clear that the soil moisture range is much better represented by the H-TESSEL scheme. This results also in a slightly better match for evaporation during the dry season. TESSEL retains too much water for this desert savannah site.
2) BERMS boreal forest
The Boreal Ecosystem Research and Monitoring Sites (BERMS) site located in a Canadian boreal Old Aspen forest (central Saskatchewan) has a high soil water retention. The data have been used for the validation of fluxes from ERA-40 in Betts et al. (2006). BERMS data represent one of the longest comprehensive time series available for land surface model verification. The data have a marked interannual variability, allowing for the evaluation of multiple time scales. Soil moisture observations are collected daily at several depths and here evaluated over the first 90 cm, which are compared to the top 100 cm of model obtained, adding the soil moisture of the first three layers. Compared to the SEBEX site, TESSEL soil texture (medium) properties are closer to H-TESSEL (medium fine, according to FAO). Still, both the absolute values and the interannual variability of soil moisture are better captured by the H-TESSEL, as shown in Fig. 8.
This is explained by a more appropriate soil texture and by the increased memory in the van Genuchten hydrology scheme. In fact, according to Fig. 2, the conductivity is greatly reduced for all the soil texture classes except coarse textured soils, leading to a much longer recharge/discharge period.
4. Offline regional and global simulations
The verification of the models’ hydrology for large domains is a complex task. This is due to both the lack of direct observations and to a composite effect of shortcomings in land surface parameterizations, which produce errors not easily traced to a single process. Global atmospheric reanalyses [ECMWF reanalysis series ERA-15, ERA-40, ERA-Interim; National Centers for Environmental Predictions–Department of Energy (NCEP–DOE) reanalysis; Japanese 25-Year Reanalysis (JRA-25)] offer an estimate of water budgets on the global domain. However, those estimates are heavily model-dependent and known to have deficiencies, either in the land surface scheme or in the surface fluxes. This limits the validity for quantitative estimates. Seneviratne et al. (2004) and Hirschi et al. (2006a) bypassed the strong dependence on the land surface model formulation used in the reanalysis system by considering the water budget of large basins using both atmospheric moisture convergence fields derived from reanalysis and observed river discharge. Water balance closure can thus be calculated without the use of atmospheric model precipitation and evaporation from the land surface model. A BSWB dataset has been gathered for major world river catchments. Considering hydrology over large domains allows for the verification of subgrid-scale parameterizations (e.g., runoff) that cannot be evaluated at single instrumented sites and for the assessment of the overall behavior of the scheme, which may be hidden in field site experiments (lack of representativity). We focus essentially on runoff and terrestrial water storage changes to validate the performance of H-TESSEL and TESSEL. Monthly time scales are considered for the largest catchments (mainly in the Northern Hemisphere). Runoff at daily time scales is evaluated on a well-known catchment experiment, Rhone-AGG (Boone et al. 2004).
a. The GSWP-2 experiment
The GSWP-2 experiment (Dirmeyer et al. 1999, 2002; Gao et al. 2004) provides a set of near-surface forcing to drive land surface schemes offline. GSWP-2 covers the period from July 1982 to December 1995. The atmospheric forcing data are provided at a resolution of 1° globally on a domain of 360 × 150 grid points and does not consider latitudes south of 60°S. The GSWP-2 data were originally based on atmospheric reanalyses (NCEP–DOE reanalysis) at 3-h intervals, corrected using observational data as detailed in Dirmeyer et al. (2002). For the current experiments, we have used the latest release of GSWP-2 atmospheric forcing based on ERA-40, where only precipitation is corrected using the Global Precipitation Climatology Project (GPCP) dataset. These fields are labeled as ERAGSWP. Temperature and pressure fields are rescaled according to elevation differences between the reanalysis model topography and that used in GSWP-2, while surface radiation and wind fields are unchanged. The state variables of surface pressure, air temperature, and specific humidity at 2 m as well as wind at 10 m are provided as instantaneous values. The surface radiation and precipitation flux represent 3-h averages, and they are kept constant over a 3-h period to ensure conservation. The instantaneous forcings are linearly interpolated in time. The integration time step for the simulation is 30 min. The 10-yr runs are aggregated into monthly climatologies and considered over large catchments.
1) Terrestrial water storage
The variation of the TWS can be expressed as
where P, E, and R are the precipitation, evapotranspiration and runoff, respectively. TWS accounts from both snowpack and soil moisture variations. As previously mentioned, Eq. (12) can be combined with the atmospheric water balance to eliminate the P and E terms, which in atmospheric reanalysis are only indirectly constrained by observations and strongly rely on model estimates. The expression obtained is
where the terms dW/dt (atmospheric total column water variations) and −∇Q (atmospheric moisture convergence) are taken from the ERA-40 reanalysis and are thus directly constrained by assimilated observations (e.g., radiosondes). The runoff R is obtained by the observed river discharge from the GRDC divided by the basin area.
The TWS estimates in the BSWB dataset, obtained from Eq. (13), are shown to correlate well with land surface soil moisture observations averaged on large domains (Seneviratne et al. 2004). Moreover, mass variations detected from the Gravity Recovery and Climate Experiment (GRACE) satellite mission correlate reasonably well with TWS (Hirschi et al. 2006b). In the GSWP-2 simulations, the TWS variations are computed directly from the monthly means of total column soil moisture variations and the snow water equivalent variations. A number of central European catchments are selected to show the effect of the revised hydrology. An improved description of the TWS change in H-TESSEL-offline can be seen in Fig. 9, both when compared to ERA-40 (in which TESSEL is used in a coupled mode) and TESSEL-offline (driven by ERAGSWP forcing). The reason for this improvement with respect to ERA-40 TWS is discussed in section 5c. The slightly larger amplitude of TWS annual cycle in H-TESSEL is mostly associated with the increased water holding capacity, as reported in Table 2.
2) Monthly runoff
The Global Runoff Data Centre operates under the auspice of the World Meteorological Organization and provides data for verification of atmospheric and hydrologic models. The GRDC database is updated continuously and contains daily and monthly discharge data information for more than 3000 hydrologic stations in river basins located in 143 countries. Over the GSWP-2 period, the runoff data for 1352 discharge gauging stations were available. A 0.5° regular-grid mean annual runoff product is provided by the GRDC (available online at http://www.grdc.sr.unh.edu/). This product is produced with an underlying precipitation climatology (Fekete et al. 2000), which is readjusted whenever the observed runoff exceeds the input precipitation. Given the large number of assumptions involved in the production of a global runoff climatology and particularly for precipitation in ungauged basins, only a qualitative comparison is appropriate. Figure 10 shows the comparison of H-TESSEL and TESSEL with the GRDC product. Improvements are mainly visible in the Northern Hemisphere, where the nonzero runoff area is significantly increased in H-TESSEL, consistent with GRDC. The limitations of the ERA–GPCP precipitation dataset (used to force the stand-alone simulations) are clearly visible in the tropical belt. In fact, a strong underestimation of rainfall [as reported in Betts et al. (2005) for the Amazon basin] produces noticeable effects on the simulated annual runoff: both H-TESSEL and TESSEL greatly underestimate the GRDC runoff rates in this areas.
A basin-scale evaluation is considered as well. It is appropriate for large basins, which are reasonably well gauged and where the time delay of the water path (routing) does not limit the validity of the comparison (e.g., excluding therefore the Amazon, where the routing plays a major role in timing the runoff). Figure 11 shows the errors of the monthly runoff evaluated on 41 different basins, listed in Table 4. RMSE and BIAS are calculated on the (TESSEL and H-TESSEL) monthly runoff against the GRDC monthly runoff estimates.
The H-TESSEL simulation shows an overall reduction of the mean BIAS for most basins (with the Ohio River basin, No. 14 in Fig. 11, being a noticeable exception). For RMSE, roughly a third of the basins (Nos. 1–14) register a net deterioration, while the majority shows a marked improvement. This behavior can be explained from the effect of snow ablation. The TESSEL snow treatment (not modified in the H-TESSEL scheme) suffers from the lack of a refreezing mechanism for water in the snowpack, which activates the runoff too quickly and produces a pronounced peak.
Although the RMSE of the monthly runoff has deteriorated for snow-dominated basins, the BIAS is generally reduced. The magnitude of the spring runoff peak is also in better agreement with observations for H-TESSEL, as shown in Fig. 12, for the Yenisei Siberian basin (No. 9 in Fig. 11). The correct timing will only be achieved by a revision of the snow scheme. The link between snow and soil moisture errors is discussed further in section 5c.
b. The Rhone-AGG experiment
Earlier evaluations of predecessors of H-TESSEL were carried out in the context of the so-called Rhone-AGG experiment (Boone et al. 2004). This experiment was designed to test the effects of spatial aggregation between finescale and coarse-scale grids for hydrological simulations of a complex terrain like the Rhône catchment valley (95 000 km2). Similar to the GSWP-2 set-up, a 3-hourly forcing was provided for the domain at different resolutions (8 km and 1°) for a 4-yr period (from 1 August 1985 to 1 August 1989). The Rhone-AGG 8-km forcing used for this experiment is based on a dense observation network, as described in Boone et al. (2004). The first year was considered to be a spin-up year. Verifying daily discharge data were also made available for a number of subcatchments and the major Rhône branch at Viviers.
At daily time scales, river routing of the modeled runoff generation greatly affects the results for catchments of the size of the Rhône basin. The RMSE and Nash–Sutcliffe scores of daily discharge data for this limited period are worse for H-TESSEL than for TESSEL (Table 5). Scores are calculated by direct comparison of modeled surface runoff generation to daily discharge measured at Viviers, and therefore are clearly penalizing the scheme with greater variability but still lacking a routing scheme. A detailed routing scheme was not available, but a crude estimate of the gridpoint-dependent delay between the modeled runoff generation and the outlet in Viviers was made by assuming a fixed channel propagation speed (100 km day−1) over the (shortest) distance between the grid point and the outlet position. Delay times were rounded to whole days (up to six days for the upper north part of the basin). This procedure (denoted by delayed output in Table 5) clearly improves the RMSE and Nash–Sutcliffe scores for H-TESSEL and deteriorates the TESSEL output.
Apart from the objective statistical scores, a clear effect of the subgrid runoff scheme on the time series of the simulated discharge is visible in Fig. 13. The high-frequency variability in TESSEL is significantly lower than in H-TESSEL (see also the spectrum shown in Fig. 14). The runoff delay from the grid points far from the outlet has a strong effect on the time series at the outlet, both for TESSEL (Fig. 13) and H-TESSEL (Fig. 14). A strong smoothing is caused by the delayed output produced by TESSEL. Apparently, peaks generated by saturated grid boxes receiving additional precipitation or melt events are compensated by grid boxes with low runoff generation, and virtually all major peaks are removed. In H-TESSEL, the delayed output causes a much better resemblance of the modeled power spectrum to the observations (Fig. 14). Without the output delay, the power at high frequencies is too high in H-TESSEL, and this overestimation is effectively removed. The contribution of high frequencies to the overall signal variance in TESSEL is clearly underestimated. Monthly discharge for the Rhône basin simulated in the GSWP-2 context resulted in a clear improvement of the RMSE score (from 0.45 mm day−1 for TESSEL to 0.26 mm day−1 for H-TESSEL) and had virtually no effect on the bias.
5. Atmospheric-coupled simulations
Global atmospheric coupled experiments are used to evaluate the land–atmosphere feedback, especially on near-surface atmospheric quantities. Various ensemble sets of multimonth integrations are performed. We refer to these integrations as climate simulations. Short-term predictions (12 h) embedded in long data assimilation cycles (covering seven months from April to October) are also used to evaluate the analysis increments, as a measure of model improvements.
a. Climate simulations
Two sets of climate simulations have been performed, using an atmospheric model configuration with 91 vertical layers and a gridpoint resolution of about 120 km (truncation T159 of a Gaussian reduced grid). A 13-month four-member ensemble for the period from 1 August 2000 to 31 August 2001 and a multiyear 10-member ensemble of 7-month runs (starting 1 April each year from 1991 to 1999) were executed. The initial conditions are provided by the ERA-40 reanalysis, and the first month is not considered for evaluation. The experiments showed a very small atmospheric response to the land surface scheme modifications. A slightly positive effect is seen in temperature at 2 m compared with ERA-40 climatology, as shown in Fig. 15. The 2-m temperature RMSE error is reduced from 1.19 to 1.08 K, while the 2-m dewpoint temperature error is unchanged. The error patterns are similar for both 2-m temperatures, and the error reduction in the HTESSEL simulation (Figs. 15b and 15d) is mostly concentrated in Northern Hemisphere continental areas.
b. Extended data assimilation experiments
Considering the land surface water budget during a data assimilation cycle, Eq. (12) becomes
where δA represents the analysis increments (in snow δSn and soil moisture δθ) added for each cycle of the data assimilation system. It is assumed that a better land surface hydrology will lead to smaller systematic increments δA. We compared the increments for two data assimilation cycles with TESSEL and H-TESSEL as land surface schemes. For H-TESSEL, the optimum interpolation (OI) soil moisture analysis (Mahfouf 1991; Douville et al. 2000; Mahfouf et al. 2000) has also been revised for consistency. The OI coefficients have been rescaled according to the ratio of water holding capacity in H-TESSEL, a function of local soil texture, to the constant value in TESSEL (see Table 2).
To focus on time scales relevant for the land surface processes, we considered a 7-month period covering the boreal summer (1 April 2006 to 31 October 2006). The analysis makes use of short-term forecast errors in 2-m temperature and relative humidity to correct soil moisture errors via a set of (OI) coefficients [see Douville et al. (2000) for the detailed formulation]. The absolute difference (H-TESSEL − TESSEL) of the soil moisture analysis increments over summer [June–August (JJA)] is shown in Fig. 16.
Generally, a reduction of the mean daily soil moisture increments is observed at midlatitudes, particularly over Europe and the central United States, where the Synoptic Ocean Prediction (SYNOP) network is dense and the OI analysis most effective. This signal confirms the positive effect of H-TESSEL over these areas, certainly when we recall that the dynamical range of soil moisture is greatly increased. Northern regions show larger soil moisture analysis increments in H-TESSEL for the month of June, associated with land surface spinup from previous snow cover and frozen soil conditions. Snow modeling errors (in particular the anticipated snow melting) contribute to larger soil moisture errors consistently with the offline results and therefore cannot be attributed to the soil hydrology. For illustration purposes, three panels are included in Fig. 16, showing that the analysis increments are mostly reduced everywhere in August (Fig. 16c) when the snow reaches its minimum extent. The interaction of snow and soil hydrology has also been highlighted by Decharme and Douville (2007) as a crucial aspect for the correct representation of the water cycle in Northern latitude basins. This feature is well known in ERA-40 reanalysis and is briefly commented in the following subsection.
c. ERA-40 analysis increments
The long ERA-40 dataset, with a frozen configuration of the data assimilation and modeling systems, provides the opportunity to produce valuable diagnostics on data assimilation increments. It turns out that the soil moisture and snow mass errors are tightly coupled at high latitudes. In ERA-40, the snow mass is analyzed from SYNOP snow depth observations, with relaxation to climatology (12-day time scale) as a constraint in data-sparse areas. Systematic soil moisture/snow mass analysis increments indicate deficiencies in the system, which can be a result of data or model errors, or a result of a suboptimal data assimilation. Figure 17 shows the monthly-mean data assimilation increments (expressed in millimeters per day) for snow and soil moisture, calculated from ERA-40 and evaluated for the decade 1986–95 (consistent with the GSWP-2 period). Results are presented in the form of a Hovmöller diagram (land-only zonal means as a function of latitude and month of the year). Focusing on northern latitudes, the snow analysis increments moving from 40°N in winter to 70°–80°N in June indicate a persistent positive correction, which is attributed to the early melting of snow. This pattern is mirrored in the soil moisture analysis, which removes the water supply. The snow depth increments and the associated soil moisture increments clearly suggest that the model is melting the snow too early, which is consistent with the offline simulations, as shown in Fig. 12. In the latter case, the snowmelt results in a runoff peak, which is one month too early compared to the observations.
At midlatitudes during the boreal summer months, positive soil moisture analysis increments are systematic (Fig. 17), which suggests that the land surface model has a tendency of running dry in summer. From budget diagnostics on the European basins, we know that the amplitude of the seasonal cycle of TWS change in ERA-40 is much smaller than observed (Fig. 9). This is also concluded by Drusch and Viterbo (2007) by comparing analyzed soil moisture with in situ observations. More insight can be gained by considering, for the same basin, the water budget terms of Eq. (14) in stand-alone simulations with GSWP2 forcing and in ERA-40 (Fig. 18). It is shown that much of the summer increments (up to 1 mm day−1) in ERA-40 are a response to the precipitation deficit in the short-range forecasts of ERA-40 (also confirmed by van den Hurk et al. 2008). It is interesting that both systems—that is, ERA-40 affected by biased precipitation and with large data assimilation increments and the stand-alone simulation with more realistic precipitation and no data assimilation—simulate very similar evaporation. Still, the amplitude of the seasonal cycle of TWS is smaller in the ERA-40 system than in the stand-alone simulations [consistent with the results of Ferranti and Viterbo (2006)]. Having similar evaporation with different soil moisture may sound contradictory. Apparently the lack of precipitation is overcorrected by data assimilation increments, suggesting that the data assimilation might be suboptimal.
The conclusion we draw here is that data assimilation increments are a valuable diagnostic. They obviously indicate a deficiency in the system, in the sense that somehow the model is not capable of matching the observations. However, data assimilation increments by themselves do not indicate the cause of the problems. In the case of ERA-40, it turned out to be particularly helpful to be able to do stand-alone simulations with observed precipitation. Both the precipitation bias and the suboptimal data assimilation contribute to the rather large soil moisture increments, which are affecting the evolution of the TWS cycle.
6. Summary and conclusions
A revised hydrology for the ECMWF land surface scheme (H-TESSEL) is tested and compared against observations from field site experiments. The revision is introduced to correct two main shortcomings of the land surface scheme: the absence of surface runoff and a globally uniform soil texture. A new dataset for soil type is included by assigning a hydrological class (up to six) to each grid cell. A revised infiltration scheme with a subgrid surface runoff description is also introduced and evaluated. In point comparisons with field site experiments, these modifications show a shift in the soil moisture range to give better agreement with observations. The soil physiographic parameters (wilting point and field capacity) associated with each soil texture produce a larger water holding capacity. In drylands, the shift of the soil moisture range gives slightly better evaporation. A boreal forest site is selected for its long time series of soil moisture observations. The interannual variability of root-zone soil moisture in a boreal forest site shows improvements with a satisfactory match to the 8-yr dataset. The improved match of soil moisture to observed values is a necessary step toward the improved assimilation of satellite data, such as microwave radiances. In conclusion, the proposed changes to the ECMWF scheme address known shortcomings without affecting the generally good performance of the land surface model in providing the lower boundary conditions to the atmospheric model.
A set of regional stand-alone experiments (GSWP-2, 1986–95) is used to evaluate the terrestrial water storage variations over the central European river basins in comparison with independent estimates based on atmospheric moisture convergence data and river discharge observations. The model global annual runoff map shows some small improvement when compared with the river discharge product. Quantitative evaluation of the runoff at monthly time scales shows a net improvement of runoff timing in relevant catchments. In basins dominated by snow, spring snowmelt is still too early because of errors in the snow scheme. The annual BIAS in runoff is reduced for most of the basins considered. Errors in snowmelt are combined with a more active surface runoff generation in H-TESSEL and therefore lead to increased RMSE in runoff. Daily time scale for the runoff are studied in the Rhone-AGG experiment (1985–89), where the main improvement is a spectrum of river discharge values closer to observations.
Atmospheric coupled verification is carried out using climate runs and long data assimilation experiments to evaluate the ensemble of changes in the land surface scheme and in the soil moisture analysis. At global scales, precipitation, snow, and soil water errors are tightly coupled to water budget issues. Standalone integrations allow us to separate errors in ERA-40 terrestrial water storage largely due to precipitation deficits in the 0–6-h forecast (and to a suboptimal soil moisture correction) from snow-related errors, which remain a feature of the new land surface scheme. The new soil hydrology appears to clearly benefit runoff and terrestrial water storage in snow-free areas. The operational implementation of a revised hydrology scheme is presented in the appendix. Future efforts to improve the annual terrestrial water balance will include adding a seasonal vegetation cycle and improving the snow scheme.
The authors wish to thank Adrian Tompkins, Thomas Jung, Paco Doblas-Reyes, and Antje Weisheimer for their help with climate runs, and Matthias Drusch and Jean-François Mahfouf for discussions on the OI data assimilation. Emanuel Dutra and Patricia de Rosnay were helpful with the independent evaluation of H-TESSEL and for valuable discussions. Deborah Salmond, Mats Hamrud, and Nils Wedi are thanked for their help with IFS code, and Jan Haseler and Joerg Urban for their help with the scripts. Aaron Boone and Joel Noilhan are acknowledged for providing the Rhone-AGG dataset and for valuable discussions. We thank the U.K. Centre for Ecology and Hydrology and the Fluxnet-Canada Research Network for making available the SEBEX and BERMS datasets, respectively. Freddy Nachtergaele and Attila Nemes are acknowledged for their help with soil texture data and valuable discussions on soil properties aggregation. Alan Betts acknowledges support from NSF through Grant ATM-0529797. Rob Hine and Anabel Bowen are thanked for their help in improving the figures.
Global Scale Implementation
The IFS at ECMWF includes the same land surface scheme that has to operate at various spatial resolutions, conserving the main features and, most importantly, the same water and energy budgets. Conservation of water storage when moving across resolution is a desirable property, although no scheme can easily achieve this, since averaging soil (as well as vegetation) properties lead to creation/destruction of information. As illustrated, the new parameterizations proposed for the soil hydrology consider a soil type (texture) that is used to assign hydrological properties and an orographic parameter that accounts for the terrain complexity.
The FAO/United Nations Educational, Scientific and Cultural Organization (UNESCO) Digital Soil Map of the World (DSMW) (FAO 2003) is available at the resolution of 5′ × 5′ (about 10 km). The FAO DSMW provides information on two levels of soil depth: 0–30 and 30–100 cm. The deep-soil layer (30–100 cm) is used to prescribe the soil texture on the soil column when moving across resolutions. The orographic runoff parameter, which determines the response of surface runoff to precipitation as a function of the complexity of the terrain, uses the standard deviation of orography. This quantity is already evaluated at high resolution and is used in the turbulent orographic form drag (Beljaars et al. 2004). In Fig. A1, the soil texture map and the standard deviation of orography are reported for the 10-day forecast operational resolution.
The soil texture and the standard deviation of orography are interpolated at various model resolutions. The standard deviation of orography can be interpolated bilinearly without special treatment. For the soil texture class interpolation, the choice of a dominant soil texture aggregation is adopted, since it largely preserves the overall number of grid points in each class.
Thus water storage properties, when moving across resolutions, is also preserved without the creation of information. Sand and clay percentages, often used in place of textural classes, do not conserve hydraulic properties when moving across resolutions. Table A1 shows the percentage of land points in each soil texture class at T159 (about 125 km) and T799 (about 25 km) resolution.
The adoption of a dominant soil texture permits an easy transition between TESSEL and H-TESSEL soil moisture, which can be obtained by a linear rescale procedure, taking into account field capacity and permanent wilting point (as detailed in http://www.ecmwf.int/products/changes/soil_hydrology_cy32r3/). The soil moisture rescale between TESSEL and H-TESSEL is introduced to preserve atmospheric effect in terms of evaporation.
Corresponding author address: Gianpaolo Balsamo, European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading RG2 9AX, United Kingdom. Email: email@example.com