The Water and Global Change (WATCH) project evaluation of the terrestrial water cycle involves using land surface models and general hydrological models to assess hydrologically important variables including evaporation, soil moisture, and runoff. Such models require meteorological forcing data, and this paper describes the creation of the WATCH Forcing Data for 1958–2001 based on the 40-yr ECMWF Re-Analysis (ERA-40) and for 1901–57 based on reordered reanalysis data. It also discusses and analyses model-independent estimates of reference crop evaporation. Global average annual cumulative reference crop evaporation was selected as a widely adopted measure of potential evapotranspiration. It exhibits no significant trend from 1979 to 2001 although there are significant long-term increases in global average vapor pressure deficit and concurrent significant decreases in global average net radiation and wind speed. The near-constant global average of annual reference crop evaporation in the late twentieth century masks significant decreases in some regions (e.g., the Murray–Darling basin) with significant increases in others.
As the earth’s whole climate system slowly changes there are likely to be greater and faster regional changes. Studies of the impacts of these changes on essential services such as fresh water supply are being made by many researchers (e.g., Harding et al. 2011) with the change in evaporation being a key aspect. Observations of large-scale evaporation over the last half century (the most studied period) are, however, not available. Consequently, models of evaporation are frequently used as an alternative. In such models the key factors that determine changes in evaporation are changes in meteorological factors such as radiation, wind speed, air temperature, and humidity.
Studies have analyzed pan evaporation data (Roderick and Farquhar 2002; Roderick et al. 2007) and reported changes in the external drivers on evaporation when there is no change in available water. In Australia these studies have demonstrated that large-scale change in wind speed (global stilling) is responsible for an observed drop in pan evaporation, although decreases–increases in radiation (global dimming–brightening) are perhaps responsible for changes elsewhere. Shuttleworth et al. (2009) demonstrated that it is not always possible to use pan evaporation to diagnose large-scale change in external drivers of actual evaporation. This is because some changes in the drivers of pan evaporation are caused by feedbacks in the atmospheric planetary boundary layer caused by altered actual evaporation in the area surrounding the pan. However, they also demonstrated that it is not possible to assume that changes in pan evaporation are equal and opposite to changes in surrounding actual evaporation, as suggested by Bouchet (1963), since changes in the variables controlling evaporation are a mixture of regional atmospheric feedbacks superposed on modified large-scale atmospheric circulation.
In their comprehensive review, Hobbins et al. (2008) point out that researchers interested in global evaporation need an accurate assessment of the external drivers of the evaporation process. However, because of nonlinearity in the relationships between the drivers of evaporation (particularly temperature) it is not possible to make such an assessment using daily average meteorological data. Instead, accurate assessment requires data that resolve the full diurnal cycle. This paper describes the creation of the Water and Global Change (WATCH) Forcing Data (WFD), a dataset that is available for the whole of the twentieth century and that resolves the full diurnal cycle. An analysis of changes in the external drivers of evaporation that is relevant to both researchers and water-resource engineers is also made.
The European Union WATCH project (www.eu-watch.org) seeks to assess the terrestrial water cycle in the context of global change in the twentieth and twenty-first centuries. A major component of the study is use of land surface models (LSMs) and general hydrological models (GHMs) to calculate changes in hydrologically important variables such as evaporation, soil moisture, and runoff (Haddeland et al. 2011). For both types of model, meteorological “forcing” (or “driving”) data (air temperature, rainfall/snowfall, etc.) are required at subdaily time steps for the LSMs and daily time steps for the GHMs. The 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) product, which provided the basis data used in the derivation of the WFD, was derived from successive short-term integrations of a general circulation model (GCM) that assimilated [via three-dimensional variational data assimilation (3D-Var)] various satellite data along with atmospheric soundings and land and sea surface observations (Uppala et al. 2005). The reanalysis procedure used to create ERA-40 merged global subdaily observations with a prior estimate based on short integrations of a comprehensive GCM, allowing for uncertainties in each, using a GCM configuration that was consistent, as opposed to the progressively refined and improved GCMs that are used in routine weather forecasting. As explained below, the WFD were derived from the surface variables of the ERA-40 reanalysis product for the period 1958 to 2001, but from reordered ERA-40 data for the period 1901 to 1957.
The several models involved in the WATCH project calculate hydrological variables using the WFD in different ways, but a key aspect of the models is the way in which evaporation is estimated (Haddeland et al. 2011). LSMs typically estimate actual evaporation by evaluating the energy balance at the subdaily time scale, whereas GHMs typically require estimates of daily-average “potential” evapotranspiration and then assess actual evaporation by adjusting this estimate to allow for the water availability. In this paper an assessment is made of changes in global twentieth-century potential evaporation independent of any specific LSM or GHM as estimated via the WFD themselves. Consideration is also given to regional variations in the selected large river basins shown in Fig. 1.
2. The WATCH Forcing Data
The WFD consist of subdaily, regularly (latitude–longitude) gridded, half-degree resolution, meteorological forcing data. The variables included are (i) wind speed at 10 m, (ii) air temperature at 2 m, (iii) surface pressure, (iv) specific humidity at 2 m, (v) downward longwave radiation flux, (vi) downward shortwave radiation flux, (vii) rainfall rate, and (viii) snowfall rate. These global data are stored at 67 420 points over land (excluding the Antarctic), with the land–sea mask used being that defined by the Climatic Research Unit (CRU; New et al. 1999, 2000) in netCDF format using the Assistance for Land-Surface Modelling Activities (ALMA) convention (see http://www.lmd.jussieu.fr/~polcher/ALMA/). Variables vi–viii are not readily interpolated and are stored at three-hourly time steps as in the basic ERA-40 data, but to save space variables i–v are stored at 6-hourly time steps with code provided to give variable-dependent interpolation to the three-hourly time step.
a. WATCH Forcing Data 1958–2001
Generation of the WFD for the late twentieth century described in detail by Weedon et al. (2010) adopted the procedures described by Ngo-Duc et al. (2005) and Sheffield et al. (2006), but with the changes summarized in Table 1. Processing involved bilinear interpolation of each variable from the 1° ERA-40 grid to the 0.5° CRU land–sea mask. To maintain consistency, elevation corrections were then made sequentially to the interpolated temperature, surface pressure, specific humidity, and downward longwave radiation (in that order, because elevation correction of later variables requires use of previously corrected variables).
In several respects the ERA-40 data product is superior to the earlier National Center for Atmospheric Research–National Centers for Environmental Prediction (NCAR–NCEP) reanalysis used in deriving other forcing datasets (e.g., Uppala et al. 2005), but the 2-m temperatures in ERA-40 are known to lack some climatic trends and to exhibit an overall bias (Betts and Beljaars 2003; Simmons et al. 2004; Hagemann et al. 2005) despite the assimilation of relevant surface observations. Comparison of diurnal extremes in near-surface temperature in the NCAR–NCEP, ERA-40, and the (more recent) Japanese Meteorological Agency (JMA) 25-yr reanalysis (JRA-25) data reveals problems in all 3 data products (Pitman and Perkins 2009), particularly with respect to minimum temperature. For this reason the monthly average interpolated and elevation-corrected temperatures from ERA-40 were also bias-corrected (Weedon et al. 2010). Because the CRU3 data (Brohan et al. 2006) were not available at 0.5° resolution for all the required observations during creation of the WFD, CRU TS2.1 gridded observations were used for this bias correction (New et al. 1999, 2000; Mitchell and Jones 2005).
The use of CRU observations for monthly bias correction inevitably incorporates inaccuracies related to creation of the gridded products. Nevertheless, the CRU interpolation methodology based on 1961–90 anomalies (New et al. 1999, 2000) includes allowance for the “correlation length” of the variables involved, and elevation corrections and inhomogeneities between stations have been adjusted while the variable station coverage through time and spatially is documented by New et al. (1999, 2000) and Mitchell and Jones (2005). Despite these limitations the CRU dataset has been widely used for investigating global terrestrial changes through the twentieth century (e.g., Déry and Wood 2005; Gedney et al. 2006; Dang et al. 2007; Piao et al. 2009).
The CRU temperature data used include some (albeit rare) inhomogeneities. Specifically, there were steplike offsets in the values that can span several years at particular sites and also some single month outliers (which were identified as being more than five standard deviations away from the 1958–2001 monthly mean). Prior to their use for bias correction, the inhomogeneities were removed from CRU data using the method of Österle et al. (2003) and single month outlier values were replaced with the local calendar-month mean (Weedon et al. 2010). Average monthly diurnal temperature ranges were also corrected using the CRU data (Weedon et al. 2010).
2) Corrections to variables other than precipitation
The relative humidity implied by the original ERA-40 temperature, pressure, and specific humidity was interpolated bilinearly to the half-degree grid following Cosgrove et al. (2003), and the resulting values then used with the elevation- and bias-corrected temperature and pressure to calculate specific humidity. Using this method maintains consistency between variables and also avoids supersaturation. CRU observations of vapor pressure were used to make monthly average checks of the values so derived, but they were not used for bias correction because this would have compromised consistency.
Using the ERA-40 data means that there is no global unidirectional bias in the WFD downward longwave radiation with respect to the average National Aeronautics and Space Administration (NASA) Surface Radiation Budget (SRB) product (Weedon et al. 2010). This contrasts with the work of Ngo-Duc et al. (2005) and Sheffield et al. (2006), where global unidirectional bias related to the NCAR–NCEP reanalysis necessitated correction via the SRB product. Comparison with selected FLUXNET data (Weedon et al. 2010) also showed that it was not necessary to make a monthly bias correction of the WFD downward longwave radiation using the SRB3 quality-controlled longwave (QCLW) product (see http://eosweb.larc.nasa.gov/PRODOCS/srb/table_srb.html) interpolated to half degree.
Downward shortwave radiation was adjusted at the monthly time scale using CRU cloud cover and the local linear correlation between monthly average (interpolated) ERA-40 cloud cover and downward shortwave radiation (Sheffield et al. 2006; Weedon et al. 2010). Troy and Wood (2009) compared unadjusted ERA-40 radiation fluxes with other reanalysis products and observations across northern Eurasia. ERA-40 does not include adjustments for the effects of seasonal and decadal variations in atmospheric aerosol loading on downward shortwave radiation fluxes (Uppala et al. 2005), although long-term changes in aerosol loading can significantly influence downward shortwave radiation fluxes (e.g., Wild et al. 2008). A correction was therefore made for the effects of tropospheric and stratospheric aerosols on downward surface fluxes of shortwave radiation using twentieth-century aerosol optical depths (AODs) taken from a GCM combined with lookup tables of radiative transfer calculations.
Distributions of tropospheric AOD at 0.55 μm for the twentieth century were taken from simulations with the atmospheric component of the Hadley Centre Global Environmental Model version 2 (HadGEM2-A; Martin et al. 2006; Collins et al. 2008). HadGEM2-A includes representation of the following tropospheric aerosol species: sulfate, mineral dust, sea salt, black carbon from fossil fuel and from biomass burning, and secondary organic aerosols (Bellouin et al. 2007). Stratospheric aerosols from volcanic eruptions were available as zonal means (Sato et al. 1993, dataset updated in 2002). Aerosol radiative effects are represented in both the clear-sky (cloud-free) portion of each GCM grid box and the portion that is cloudy. Thus, the calculations made on the GCM grid and interpolated to half degree provided correction to clear-sky downward radiation that accounted for the direct and indirect effect of aerosols in the troposphere and the direct effect in the stratosphere, and also for the effect of aerosols on cloudy-sky downward radiation in the troposphere (Weedon et al. 2010). These corrections assume that stratospheric aerosols do not interact with tropospheric clouds to influence cloudy-sky radiation fluxes, and there is also no allowance for indirect effects of aerosols on ice clouds (cirrus) in the stratosphere. The aerosol load-corrected shortwave radiation was compared to the SRB version 3 quality-controlled shortwave (QCSW) product and both datasets were validated against FLUXNET observations; the comparison showed (Weedon et al. 2010) that it was not necessary to bias-correct the WFD downward shortwave radiation using the SRB3 QCSW product.
3) Corrections for rainfall and snowfall
The generation of the precipitation data for the WFD involved six steps (Weedon et al. 2010): 1) bilinear interpolation, 2) combining rainfall and snowfall totals while retaining the rainfall/snowfall ratio for each location and time step, 3) adjusting the number of “wet” (i.e., rain or snow) days per month to match the CRU TS2.1 observations, 4) adjusting the monthly precipitation totals to match version four of the Global Precipitation Climatology Centre full product (GPCCv4), 5) reassigning the precipitation into rain and snow using the original ratio, and 6) adjusting the monthly totals using gridded average precipitation gauge corrections (separately for rainfall and snowfall).
The GPCCv4 full data product used in step 4 is based on gridded precipitation gauge measurements comparable to the CRU totals (i.e., they exclude satellite information and do not include gauge corrections; T. Fuchs 2008, personal communication). This observational dataset was chosen for adjusting monthly precipitation totals rather than CRU TS2.1 totals because their station coverage is much better, particularly at high latitudes and for the end of the twentieth century (Rudolf and Schneider 2005; Schneider et al. 2008; Fuchs et al. 2009; see http://gpcc.dwd.de/). Exploratory precipitation processing using the CRU totals for correction instead of GPCCv4 had revealed minor differences during the boreal winter (December–February) and major differences in northeast India–Bangladesh and northern Amazonia during boreal summer (June–August; Weedon et al. 2010).
The method adopted for wet-day correction is the main difference in the derivation of previous precipitation forcing datasets. Ngo-Duc et al. (2005), for example, did not correct wet days, whereas Sheffield et al. (2006) used a statistical correction (Sheffield et al. 2004) that was designed to cope with spurious standing wavelike patterns in the high northern-latitude wet-day characteristics of the NCAR–NCEP data. However, the Sheffield et al. correction meant that spatial continuity of individual precipitation events was sometimes compromised (see Fig. 7 of Sheffield et al. 2004), and it also required the adjustment of several associated variables when wet days were “created” to match the CRU data.
The main weakness with ERA-40 precipitation is the presence of too many wet days in the tropics (Betts and Beljaars 2003; Hagemann et al. 2005; Uppala et al. 2005) rather than spurious standing wave patterns. The approach used to redress this weakness was to compare the number of wet days in a particular month at each half-degree grid square with the CRU data. When and where there were too many wet days in the interpolated data (specifically two days or more than the CRU count), the number of days with precipitation in the month was reduced by progressively setting the rainfall–snowfall rate to zero on the day with the lowest daily total precipitation until the number of wet days matched the CRU count. Resetting of the precipitation rate was made without reference to the associated specific humidity.
This method for wet-day correction has the advantage that, because only the smallest daily totals are reset, the spatial continuity and coherence of significant (nondrizzle) frontal precipitation across grid boxes is not compromised. This is important in the context of the WATCH project because it means that large-scale (multigrid box) hydrological modeling remains meaningful at the daily scale. For locations where there were too few wet days per month relative to the CRU observations, no changes were made, thus avoiding the need to artificially modify downward shortwave, specific humidity, and 2-m temperature on dry days to make them consistent with conversion to wet days (cf. Sheffield et al. 2006).
The correction method just described was successful in that the number of tropical wet days was adjusted to match the CRU data and the adjustment of precipitation totals based on GPCCv4 totals is not problematic. However, for the (very few) locations and times when there were too few wet days in the interpolated ERA-40 data, the adjustment of monthly precipitation totals sometimes implied extraordinarily high precipitation rates, and it was expedient to limit these “outlier” rates to a rate corresponding to the 99.999% lognormal probability precipitation rate for the relevant calendar month and grid box (Weedon et al. 2010). As a result, some precipitation totals are less than the GPCCv4 totals in the WFD in a few locations and months. In a small number of grid boxes and some months precipitation rates are close to zero in the 1958–2001 ERA-40 data. The monthly bias correction then had the effect of increasing these rates such as to imply there was spurious background drizzle between more normal precipitation events. In semiarid areas this is inconsistent with local climatic conditions but, fortunately from the point of view of hydrological modeling, this spurious low-level background precipitation is not significant.
Once the number of wet days and precipitation totals had been adjusted, the rainfall and snowfall proportion at each time step and grid box were assigned to the ratio of rain and snow originally diagnosed by the ERA-40 reanalysis (i.e., step 5). This means that the full atmospheric profile is involved is allocating precipitation to rain and snow rather than (say) simply using a threshold of 0°C in 2-m temperature. The subsequent precipitation gauge catch correction used separate average calendar monthly catch ratios for rainfall and snowfall rates at each half-degree grid box taken from Adam and Lettenmaier (2003), who originally provided either rainfall or snowfall catch ratios for each calendar month and grid box. No attempt was made to adjust precipitation rates to allow for the effects of orography (cf. Adam et al. 2006).
Part of the validation process for the WFD involved use of FLUXNET data (www.fluxnet.ornl.gov/fluxnet/), which were obtained (with permission) and then gap-filled for selected years at seven sites (see Fig. 1) (Persson et al. 2000; Aubinet et al. 2001; Araújo et al. 2002; Suni et al. 2003; Meyers and Hollinger 2004; Grünwald and Bernhofer 2007; Urbanski et al. 2007; Göckede et al. 2008). This selection of sites allowed direct comparison of data from the mid-1990s to 2001 (consequently restricting the geographic availability of data principally to Europe and North America) and included a variety of latitudes and climatic regimes and a variety of land-cover types and elevations.
Weedon et al. (2010) illustrate time series of several variables and also provide spatial comparisons of (a) the seasonal averages of the vapor pressure implied in WFD with data from CRU, (b) WFD downward longwave and shortwave fluxes with bias-corrected versions using SRB satellite averages, and (c) WFD precipitation with a bias-corrected version that used CRU monthly totals rather than the GPCCv4 monthly totals. The validation studies discussed here are restricted to consideration of snow–rain transitions, statistical comparison of time series, and illustration of the time series of temperature and precipitation.
The subsidiary figures in Fig. 2 compare the proportion of snowfall relative to total precipitation as a function of near-surface temperature for flux tower sites (excluding snow-free Manaus) with the corresponding proportion at equivalent half-degree grid squares in the WFD. These figures illustrate data only when precipitation rate (snowfall plus rainfall) exceeds 0.5 mm h−1; consequently, a snowfall–precipitation ratio of 0 indicates precipitation that is exclusively rainfall rather than 0 precipitation. When flux tower observers arbitrarily assigned the proportion of snow to be exactly ⅓, ½, or ⅔ of the total precipitation, these ratios were not deemed reliable and were excluded from Fig. 2.
Figure 2 shows that in both the WFD and the (original and 3-h aggregated) flux tower observations, the transition between snow and rain is not well defined by using a 0°C threshold in 2-m temperature (shown as vertical gray lines). In the flux tower observations rain alone (snow/precipitation = 0.0, precipitation ≥ 0.5 mm h−1) often occurs below this threshold, while snow alone (snow/precipitation = 1.0) also occurs above this threshold. Interestingly, between −15° and −2°C the WFD (and ERA-40 reanalysis) rarely has precipitation that is exclusively rainfall or snowfall, and in the original flux tower data a mixture of rain and snow is also fairly common. The proportion of half-hourly flux tower data that imply mixed rain and snow depends on latitude. At Hyytiala (61.85°N) 16.6% of the data are mixed phase precipitation whereas at Bondville (40.0°N) just 1.9% are mixed phase, although these percentages should be considered minima because the artificially defined sleet–wet snow observations (ratios of exactly 0.5, 0.333, and 0.666) were excluded from the figure. Overall the results indicate that using the proportions of rain and snow indicated by the WFD in hydrological modeling is likely to be more reliable than assigning a water phase based on a 2-m threshold temperature (cf. Table 1 in Haddeland et al. 2011).
Table 2 gives the squared correlation coefficient (r2, which indicates the proportion of variance shared by the two time series), the root-mean-square error (RMSE), the mean bias error (MBE, i.e., the mean data point differences), and the lag-1 autocorrelation (ρ1, the 1 time step serial dependence) between 3-h FLUXNET data and the WFD. The lag-1 autocorrelation characterizes the “red” noise (nonregular) component of time series—smoothly varying data have a value of ρ1 near 1.0 whereas very noisy/erratic data have a value near 0.0. This parameter was determined using the robust spectral-fitting method of Mann and Lees (1996) because large-amplitude regular components such as diurnal and annual cycles can cause a positive bias. Correlation coefficients were calculated having removed the lag-1 autocorrelation, which otherwise positively biases the calculation, via prewhitening of the time series (i.e., Xt,pw = Xt − ρ1Xt−1, where Xt,pw represents the prewhitened value of the time series at time t; e.g., Ebisuzaki 1997). For precipitation and shortwave radiation the number of data points used in the calculation of Student’s t, used to assess the significance of the correlations, was reduced by excluding from consideration times of zero precipitation and night time values, respectively.
It should be recognized that data in the WFD represent half-degree grid box area averages but FLUXNET data represent very much smaller sensor “footprints” (Göckede et al. 2008). The correlations between these two sources of data are highly significant for all locations and variables, with the notable exception of precipitation at Manaus and Harvard Forest, largely because of the very large sample sizes (Table 2). However, several variables sometimes have large shared variance, specifically 2-m temperature (r2 = 0.21–0.64), surface pressure (r2 = 0.09–0.37), downward longwave radiation (r2 = 0.05–0.48), and downward shortwave radiation (r2 = 0.65–0.84). Conversely, correlation of prewhitened specific humidity is low at all sites (r2 = 0.03–0.12) although RMSE and mean bias errors are low compared to the means.
In Fig. 3 daily average WFD 2-m temperature is overlaid (in gray) on half-hourly flux tower values (in black). The daily 2-m temperature tracks the center of the half-hourly (hourly for Harvard Forest and Manaus) field data well, indicating that the WFD capture local daily-to-monthly (synoptic) meteorological variability as well as the seasonal cycles. The general similarity in values at the different spatial scales of the WFD and the field observations is symptomatic of the long spatial correlation length of temperature (New et al. 2000).
The only selected flux tower that is located in an area of predominantly convective rainfall is at Manaus in Amazonia. Although the number of wet days each month and monthly total precipitation had been adjusted in the WFD, at the three-hourly time scale the development of cloud and the occurrence of convective rainfall in the reanalysis for this site only poorly match the flux tower observations, even when the latter are aggregated to give three-hourly values. At the other flux tower sites considered, rainfall and snowfall associated with frontal systems in the reanalysis are more likely to match field observations at the daily to monthly time scales because the probability of precipitation is partly influenced by assimilated observations (such as atmospheric pressure). Overall the correlations for precipitation are low (r2 = 0.000–0.046) and the root-mean-square error is large. Mean bias error indicates overall mismatch in values over the full duration of the data in Table 2, and the assertion that the match is better at longer time scales is supported by the low absolute values of the MBE compared to the mean precipitation at all locations. Additionally, Fig. 4 shows that at several flux tower sites both the occurrence and intensity of daily precipitation in the WFD show a good match to daily average observations (e.g., Hyytiala and Harvard Forest).
The r2 of the prewhitened time series is below, and sometimes far below, 0.1 for wind speed at all sites except Vielsalm, and at Bondville the mean bias error for wind speed is especially large compared to the mean. This is likely to be because the Bondville flux tower is located in an area of crops while the reanalysis treats the full grid square as being forest. As a result, generally high and very variable observed winds are being inappropriately compared with generally low and much less variable modeled forest-cover winds.
At Collelongo correlations are low in comparison with other sites for 2-m temperature, specific humidity, and downward longwave radiation, and the mean bias error is also high for these variables. A likely contribution to these discrepancies is that the flux tower site is 564 m higher than the grid box average elevation (Table 2). This affects the 2-m temperature (via the environmental lapse rate) and also surface pressure, and these two variables in turn influence specific humidity and downward longwave radiation and hence the mean bias error. It is likely that local topographic factors also led to a mismatch (i.e., low correlations) between the flux tower weather and the grid-square average reanalysis results.
The RMSE for downward shortwave radiation is fairly high (~90 W m−2) at all sites and especially so at Manaus (109 W m−2). This is expected because convective clouds are difficult to model correctly in GCMs so there is likely to be a large mismatch with the field observations at the 3-h scale. However, absolute mean bias errors are acceptable (2–23 W m−2; Table 2) and the correlations are high since the CRU fractional cloud cover was used to correct mean downward shortwave radiation in the WFD at the monthly scale [see section 2a(2)].
The lag-1 autocorrelations show an impressive level of agreement at all localities for all variables with the exception of wind speed and precipitation. The reanalysis wind speed often has a higher lag-1 autocorrelation than observations (i.e., the variability between the three hourly time steps is too low), although for some unknown reason the opposite is true at Hyytiala. At all the sites the precipitation lag-1 autocorrelation is always very much higher in the WFD than for observations, indicating that, compared with reality, there is too much serial dependence (“memory” or “inertia”) in the generation of precipitation in the GCM, at least at these sites.
b. WATCH Forcing Data 1901–57
To allow modeling of hydrological processes in the WATCH project for the full twentieth century, forcing data are required for 1901–57, but prior to 1958 reanalysis data from ERA-40 are not available. It is therefore necessary to create a data series of key variables for each grid box that have appropriate characteristics in terms of their diurnal to monthly variations. These data were generated using reordered ERA-40 data a year at a time rather than by using a “weather generator.” This approach has the advantage that it ensures spatial coherence of frontal rainfall and snowfall events across grid boxes, which is very important for hydrological modeling of large river basins but which is difficult to ensure in data created using a weather generator. Additionally, the procedures adopted guarantee that the ensuing data has the same temporal variability (diurnal, submonthly variations), the same autocorrelation characteristics (serial dependence from subdiurnal to yearly scales), and the same covariance relationships between variables as during the ERA-40 interval. The procedures used to create the WFD for the period 1901–57 are described below.
1) ERA-40 data assignment
Separate years of ERA-40 data were extracted in their entirety to provide the basic data. The extraction order used (see Table A1) was random, based on the ran1 algorithm of Press et al. (1992), subject to the following constraints:
Years of ERA-40 data were extracted in random order and assigned in random order without replacement to the years 1901–57 until all 44 of the ERA-40 years from 1958–2001 had been extracted.
The 13 remaining years of required data needed were assigned again in random order without replacement until all 57 years had been allocated ERA-40 data.
In the selection process only leap years were assigned to leap years and only nonleap years were assigned to nonleap years.
This selection procedure ensures that as a global average, the statistical characteristics (e.g., overall frequency of daily to seasonal extremes) of the assigned data for 1901–57 are the same as for 1957–2001. Note that the timing of particular weather events (e.g., exceptional precipitation) is certainly not correct at any particular site, as would also have been the case had a weather generator been used.
2) Data adjustments
Exactly the same initial processing steps were applied to the 1901–57 basic data as to the 1958–2001 data (i.e., bilinear interpolation and sequential elevation corrections). The same adjustments of monthly averages (i.e., including the corrections for discontinuities and outliers and diurnal temperature range in the CRU data) were applied to 2-m temperature prior to an elevation correction of surface pressure, specific humidity, and downward longwave radiation. Downward shortwave radiation was again adjusted using the CRU cloud-cover observations, and the effects of seasonal and long-term atmospheric aerosol loading on downward shortwave radiation appropriate for 1901–57 were applied. Total precipitation was also again adjusted using the 1901–57 CRU wet days and the GPCCv4 product monthly precipitation totals prior to making separate rainfall and snowfall gauge-catch corrections.
An important factor to consider in the use of monthly bias correction of the pre-1958 data is the variable temporal and spatial coverage of the CRU and GPCCv4 meteorological station network. This has been documented by New et al. (1999, 2000), Mitchell and Jones (2005), and Fuchs et al. (2009; see http://gpcc.dwd.de/). In general the station coverage is worst prior to 1950 especially for precipitation gauges and cloud-cover observations. The regions with the most limited station coverage prior to 1950 are northern central South America, southwest China, the Sahara and central Africa, the Saudi peninsula, and high northern latitudes in Canada and Russia. For specific months and variables, for those grid boxes that have too few meteorological observations for reasonable interpolation, CRU substitutes the local monthly 1961–90 climatological average.
3) Removal of year-end discontinuities
At each grid box, reordering of complete years of ERA-40 data frequently led to year-end discontinuities in wind speed, 2-m temperature, surface pressure, specific humidity, and downward longwave radiation. This was mitigated by applying a “ramp” in the average daily values for these variables between 1 and 5 January for each year from 1902 to 1957. The mean daily values of variables at each grid box were found for 6 January and for 31 December of the preceding year (values on these days were left unchanged). Based on these, ramp adjustments were applied so that the moving-window, mean 24-h values for 1–5 January changed linearly at each 3-h time step. In this way the mean weather in one year adjusted to the mean weather in the next year over a 5-day period, this period being chosen to approximately correspond to the typical transit time of frontal systems, and so that introducing the ramp does not greatly bias the monthly average weather in January. Similar ramps were applied to the last 5 days of December 1957 data to allow a smooth transition between the pre-1958 and the original ERA-40-based 1958–2001 data.
In the case of 2-m temperature, the monthly adjustments to the CRU average temperature and diurnal temperature range were reapplied after creation of year-end ramps so that the ramped temperature agreed with the January CRU monthly averages. No year-end ramps were applied to the rainfall, snowfall, or downward shortwave data because these variables change greatly from day to day largely in response to cloud cover, and imposing a ramp in the daily values for these variables is therefore unrealistic.
3. Estimation of reference crop evaporation
To estimate actual evaporation, GHMs typically first calculate an estimate of potential evapotranspiration (PET), which is often based on either the Penman–Monteith equation (Monteith 1965) or the Priestley–Taylor equation (Priestley and Taylor 1972). This calculation seeks to characterize the evaporation (or latent heat) that might be expected from a hypothetical well-watered vegetation–soil surface that is subject to the ambient meteorological forcing variables. Models then estimate the actual evaporation as a proportion of the PET based on the land cover present and the availability of moisture in the soil or on the canopy. Thus PET can always be estimated, even for hot and cold deserts where there is little chance of significant actual evaporation because there is limited moisture available.
Changes in PET implied by the WFD from 1901 to 2001 were evaluated by calculating daily average values, but 3-h time steps of the WFD were used in this calculation because net longwave radiation and saturation vapor pressure vary nonlinearly with temperature. In humid conditions the Priestley and Taylor (1972) equation is sometimes used in GHMs (e.g., Haddeland et al. 2011) to make an estimate of potential evaporation, hereafter called PETPT (in units of W m−2); thus,
where Δ is the rate of change of saturated vapor pressure with 2-m temperature, γ is the psychometric constant, and α is a factor, usually set to 1.26 (Priestley and Taylor 1972), that apportions the available energy (A) between sensible heat and latent heat from saturated land surfaces. Assuming zero net daily ground heat flux (Allen et al. 1998), at daily time scales the available energy is usually set equal to the net radiation given (Shuttleworth et al. 2009) by
where a is the albedo (often set as 0.23 for vegetated surfaces), S is the downward shortwave radiation flux, and Ln is the net longwave (upward minus downward) radiation flux.
The Penman–Monteith equation (Monteith 1965) provides an opportunity to make an estimate of potential evaporation that allows for both the influence of available energy and atmospheric humidity on evapotranspiration through vapor pressure deficit (VPD) and wind speed. For this reason it is appropriate not only in humid but also in arid and semiarid climates. Shuttleworth (2006) and Shuttleworth et al. (2009) discussed the historical basis of the Penman–Monteith equation and practicalities of its calculation. Allen et al. (1998) specified a version of the Penman–Monteith equation that is now widely adopted as providing an estimate of evaporation from a “reference crop” (i.e., from a hypothetical, well-watered, 12-cm-high grass crop) by defining specific values of the resistances that appear in the Penman–Monteith equation. Thus, to obtain estimates of reference crop evaporation rate, hereafter referred to as PETrc, the surface resistance rs is specified as being 70 s m−1 and the aerodynamic resistance ra (in s m−1) as
where u2 is the 2-m wind speed (derived from the WFD 10-m wind speed by multiplying by 0.749; Allen et al. 1998).
The vapor pressure deficit is given by
where e is the vapor pressure and esat the saturation vapor pressure. Using ra and rs specified for the reference crop, the version of the Penman–Monteith equation that provides an estimate of PETrc in W m−2 (Shuttleworth et al. 2009) takes the form
Thus the calculation of PETrc requires use of six of the eight WFD forcing variables. The reference crop is defined to be always well watered and of limited extent, so that its presence does not significantly impact the value of the grid-average forcing variables, which are in part determined by the true area-average actual evaporation rate. If actual observations are used as forcing variables, the effect of area-average evaporation is presumably reflected in their values. However, if the forcing variables are in part derived from reanalysis data, it is implicitly assumed that the model used to calculate these (ERA-40) reanalysis data correctly calculates area-average actual evaporation, and its dependence on soil moisture. This assumption may not always be true in some regions and in some atmospheric conditions. In the following, PETrc and PETPT are compared as alternative estimates of potential evapotranspiration and have been converted to equivalent depth of evaporated water (in mm) for ease of comparison with modeling results (e.g., Haddeland et al. 2011). Lu et al. (2005) investigated a selection of radiation-based or temperature-based PET methods, adopted where the full range of observed meteorological variables is not available, and rated their performance against Food and Agriculture Organization (FAO) reference crop evaporation used as a standard. Recently Kingston et al. (2009) compared a variety of methods for evaluating potential evapotranspiration globally under climate change.
4. Global reference crop evaporation
Figure 5a shows the average cumulative PETrc per year calculated from the WFD for 1979–2001. In arid areas such as the Sahara Desert, the calculated value of PETrc far exceeds the actual evaporation (Jung et al. 2010) from natural surfaces. In fact, the areas in Fig. 5a where average PETrc exceeds 1500 mm yr−1 correspond well to the hot desert areas of the globe. As mentioned earlier, PETPT is arguably an estimate of potential evaporation that is reliable in humid areas, although it has been used in this way elsewhere in GHMs (Haddeland et al. 2011). To demonstrate the discrepancy between these two alternative estimates of potential evapotranspiration, Fig. 5b shows PETPT with the same scale as Fig. 5a. This figure clearly shows that PETPT can differ locally by more than 1000 mm yr−1 and confirms the findings of Kingston et al. (2009). In part this explains why in the Water Model Intercomparison Project (WaterMIP) exercise (Haddeland et al. 2011), which used the WFD for the period 1985–99, the GHMs using PETPT (participating alongside LSMs and GHMs using PETrc) contributed to the wide scatter in the model results for arid areas such as the upper Niger River basin, the Orange River basin, and the Murray–Darling River basin (see Fig. 6 of Haddeland et al. 2011).
Figure 6 shows changes in the global, area-weighted, annual average, cumulative PETrc during the twentieth century derived from the WFD. The gray zone around the average values indicates the 95% confidence interval of the mean assessed across all grid boxes. This uncertainty does not include assessment of the uncertainties due to the generation of the gridded CRU data for monthly bias correction. Table 3 documents the linear trends in PETrc and associated variables and their significance as assessed from the distribution of mean values around the regression, but not their uncertainty due to uncertainties in the CRU data. Trends over the period 1901–57 are calculated separately from those over the period 1958–2001. This is because, by randomizing the order of the ERA-40 basis data, the process used to create the WFD before 1958 removes the interannual dependency of variables that were not subsequently bias-corrected (wind speed, surface pressure, specific humidity, and downward longwave radiation). This change in character of interannual variations in the WFD prior to December 1958 is reflected in the more erratic changes in PETrc and other variables in Fig. 6 relative to the more smoothly varying changes after January 1958. It is also reflected in the fact that the lag-1 autocorrelation of global annual PETrc is 0.30 before 1957 and 0.64 afterward.
Throughout this paper linear trend significance is assessed using a Student’s t test in which the lag-1 autocorrelation is used to estimate the (lower) effective number of independent data points in order to allow for the influence of the serial dependence of the time series (Zwiers and von Storch 1995; von Storch and Zwiers 1999). Based on these criteria the trend in global annual PETrc from 1958–2001, which is −0.51 (±0.20) mm yr−1 per year, is statistically significant (Table 3). However, there is no significant trend in global PETrc calculated from the WFD from 1901 to 1957.
The lack of trend globally in the earlier part of the century could be a genuine phenomenon or it may in part reflect the procedure used to generate these data by use of randomized individual years of ERA-40 basis data. Although there are increases in 2-m temperature incorporated into the WFD (1901–57) via bias correction, it is likely that the lack of monthly bias correction of wind speed, surface pressure, specific humidity, and longwave radiation meant that PETrc does not incorporate climate change trends due to the randomization of the individual years of ERA-40 basis data. Potentially the use of future early twentieth-century reanalysis data could help recover possible interannual variability in PETrc. Additionally, in those locations where there were insufficient meteorological stations for interpolation prior to 1950, CRU substituted monthly 1960–91 climatology [as discussed in section 2b(2)]. In such locations the use of CRU-substituted climatological averages in bias correction, rather than real observations, will have further led to removal of any decadal and longer trends in PETrc.
The interannual variations in global PETrc are very large compared to the statistically significant linear decrease over the period 1958–2001, and they appear to have some relationship to VPD (the correlation between VPD and PETrc has r2 = 0.59, N = 44, and P < 0.001). In Fig. 6 the values of PETrc and VPD are both noticeably higher during the period 1958–73 than during the remainder of the 1958–2000 period. Uppala et al. (2005) discussed problems with the use of observations, resulting in surface pressure in the early years of the reanalysis for the periods 1958–72 and 1973–76, which they assessed as being higher and lower, respectively, than for the period 1978–2001, when use of satellite data led to more stable, better-constrained values. Because surface pressure was not bias-corrected in the WFD it is possible that the deviations in global VPD from 1958 to 1978 shown in Fig. 6 are a symptom of this feature in the ERA-40 reanalysis. There is certainly a striking similarity between features shown in our Fig. 6 and in Fig. 10 of Uppala et al. (2005). The variations in VPD in Fig. 6 are necessarily also reflected in PETrc [Eq. (5)].
There are statistically significant increases in global VPD and also statistically significant decreases global net radiation and wind speed over the period 1979–2001 (Table 3). Despite the fact that these variables have an important influence on evaporation, global PETrc shows no statistically significant change over this period. As expected, 2-m temperatures also increase substantially over 1979–2001 (see Fig. 6 and Table 3). The lack of change in PETrc over this time is presumably because of the counteracting influences of changes in other contributing variables: VPD, net radiation, and wind speed.
5. Regional reference crop evaporation
Figure 1 shows the location of eight of the large river basins that are of special interest for hydrological modeling in the WATCH project. In this study trend analyses were made for PETrc and associated variables for all eight basins (Table 4) and are illustrated for four of them (Figs. 7 and 8).
Figure 7 shows that PETrc is relatively low in the Amazon and Congo basins and also agrees fairly well, in terms of annual average, with PETPT because they lie in humid areas. In Amazonia, interannual variations in PETrc and VPD are similar to the global variations shown in Fig. 6, and PETrc had no significant trend between 1979 and 2001 although wind speed decreased significantly (Table 4a). There was also no trend in PETrc in the Congo basin from 1979 to 2001, although VPD increased significantly and there were significant decreases in wind speed and net radiation (and hence in PETPT; Table 4b).
By contrast, the Niger and Murray–Darling basins (Fig. 8) have relatively high PETrc that far exceeds PETPT because these are in arid regions. In the Niger River basin the only variable illustrated in Fig. 8 to have a significant trend over the period 1979–2001 is the wind speed (decreasing; Table 4g). In the Murray–Darling basin there is a significant decrease in PETrc over the period 1979–2001, which is probably associated with a significant decrease in VPD (Table 4d). Given the lack of a global trend in PETrc in the period 1979–2001 (see Fig. 6; Table 3), the significant downward trend in the Murray–Darling basin is clearly compensated by simultaneous increases elsewhere, providing a reminder that global changes are an amalgam of conflicting regional changes.
6. Global and regional precipitation
As previously explained, monthly precipitation totals were established for 1901–2001 by combining GPCCv4 observed totals with wet-day corrections, ERA-40 rainfall–snowfall proportions, and precipitation gauge catch corrections. Consequently, in theory interannual trends in precipitation in the WFD in the period 1901–57 are based on observations, and the reordering of ERA-40 basis precipitation data prior to bias correction does not destroy evidence of climatically significant change. However, the coverage of precipitation gauges in the early part of the century is very sparse in many regions prior to 1950. Consequently, global trends in precipitation prior to 1957 have not been assessed because the variable station coverage may have caused substantial bias in the global average values. After 1958 there is still much sparser coverage in the observing network at high latitudes (e.g., Mackenzie and Lena River basins) and in low latitudes (e.g., Amazon and Congo River basins) in comparison to midlatitudes (New et al. 1999, 2000; Fuchs et al. 2009). Consequently, the assessments of regional trends discussed below are likely to be less reliable than for midlatitude regions with the highest density of observation networks.
The WFD show no significant changes in precipitation from 1958–2001 (Fig. 9; Table 3). Globally precipitation slightly exceeds actual evaporation over land (Trenberth et al. 2007). However, PETrc is slightly larger than precipitation (shown respectively as dashed lines and lines with plus symbols in Fig. 9), although this is not surprising given that over large areas of the land surface the average reference crop evaporation exceeds average actual evaporation (cf. our Fig. 5a and Fig. 1a of Jung et al. 2010).
High-latitude cold river basins such as the Mackenzie and Lena River basins exhibit very large interannual variations in total precipitation (including snowfall), which is not seen in PETrc although the average values are similar (Fig. 9). The Mackenzie basin had decreasing snowfall from 1958–2001, but no changes in rainfall (Table 4e). The Lena basin had no significant change in either form of precipitation late in the century. The Amazon and Congo basins are relatively humid and precipitation (which is almost exclusively rainfall) far exceeds PETrc, leading to substantial runoff (Fig. 9). Rainfall apparently remained approximately constant in the Amazon basin. There were no significant changes in precipitation in the Ganges–Brahamaputra basin over 1958–2001, but there were significant decreases in the Congo River basin over the same period (though not over 1979–2001; Table 4).
Basins in arid regions such as the Murray–Darling River basin have relatively low precipitation that is almost entirely rainfall, with very large interannual variations relative to the mean (Fig. 9). These have significant implications for water resources management. In such basins potential evapotranspiration also greatly exceeds precipitation (Fig. 9; Table 4). However, the WFD show no significant trends in cumulative annual precipitation in the Murray–Darling and the Orange River basins over 1958–2001 (Table 4). Note, however, that this trend analysis does not investigate interannual changes in precipitation intensity, which may be important. The Niger River basin apparently had significantly decreasing rainfall over 1958–2001, though not over 1979–2001 (Table 4g).
This paper describes the Watch Forcing Data (WFD) created at half-degree resolution for the purpose of driving LSMs and GHMs through the twentieth century. For the period 1958–2001 the WFD can be considered to provide a good representation of real meteorological events, synoptic activity, seasonal cycles, and climate trends. The WFD for the period 1901–57 were constructed to have similar subdaily to seasonal statistical characteristics (including averages, extremes, covariance between meteorological variables, and subdaily to seasonal autocorrelation) as for the period 1958–2001. For the period 1901–57 the WFD can therefore be used to characterize early twentieth-century subdaily to seasonal hydrological statistics, but they do not represent particular historical events. There is a lack of interannual–decadal variability in PETrc for 1901–57 despite the trends in 2-m temperature introduced by bias correction as a result of 1) the randomization of the ERA-40 data used in construction and 2) lack of bias correction of wind speed, surface pressure, specific humidity, and downward longwave radiation, combined with in some regions with 3) the substitution of climatology for observations in some bias-correction data (especially cloud cover) as a result of limited observations. Potentially, effort directed toward bias correction variables (point 2) and/or use of future 1901–57 reanalysis products will alleviate these shortcomings. Nevertheless, because they are bias corrected and based directly on reanalysis, the WFD for the period 1958–2001 do include observed climatological trends in monthly to interannual changes in 2-m temperature, downward shortwave radiation, and precipitation.
When making the wet-day corrections, care was taken to avoid destroying the spatial coherence of significant precipitation events associated with frontal systems that occur across several half-degree grid squares. The WFD precipitation data also preserve the same mixture of rainfall and snowfall as in the original ERA-40 reanalysis rather than using a simplistic rain/snow threshold based on 2-m temperature. Validation against flux tower data aggregated to 3-h time steps shows that the WFD are least satisfactory in terms of describing subdaily variations in precipitation, but at monthly and longer time scales most variables show a very good level of agreement with field observations despite the difference in the spatial scales to which the WFD and flux station data relate.
Globally (excluding Antarctica), rainfall and snowfall on land remained approximately constant from 1958–2001, but is difficult to assess prior to 1957 because of inadequate and variable gauge coverage and after this time there are several areas where the changes inferred here may be biased by inadequate gauge station coverage. Snowfall apparently decreased in the Mackenzie Basin in the period 1957 to 2001. Rainfall decreased in the Congo and Niger River basins after 1958 (Table 4). There were no significant trends in precipitation in the Ganges–Brahmaputra, Orange, or Murray–Darling River basins in the twentieth century although no account was taken of interannual changes in the intensity of precipitation events in the trend analysis used.
Recognized problems with the global average surface pressure in the ERA-40 reanalysis in the period 1958–78 may well have affected calculations of global average VPD from the WFD and thence global average PETrc, and interannual variations in these two variables over this time period are probably spurious. The interannual variations in VPD and PETrc in the Amazon basin (but not the other basins studied) appear to reflect the same problems as the global data in the period 1958 to 1978.
Globally, annual average PETrc calculated using the WFD exhibits no significant change over the period 1979 to 2001 despite simultaneous significant increases in VPD and simultaneous significant decreases in net radiation and wind speed. However, the lack of overall change in global PETrc shrouds conflicting regional changes. There was, for example, a significant decrease in annual average cumulative PETrc in the Murray–Darling basin that was associated with an increase in VPD.
We thank Jan Polcher for suggesting use of the ERA-40 data as the basis for the WFD for the period 1901 to 1957, and Nigel Arnell for suggesting that corrections should be made to both rainfall and snowfall for each month. We acknowledge permission to use flux tower observations given by Steven Wofsy, Alesandro Araújo, Celso von Randow, Bart Kruijt, Christian Bernhofer, Thomas Grünwald, Timo Vesala, Giorgio Matteucci, and Marc Aubinet. This research was undertaken within the EU FP6 project WATCH (Contract 036946). GPW, NB, OB, and MB were supported by the Joint DECC/Defra Met Office Hadley Centre Climate Programme (GA01101).
ERA-40 Data Extraction Order
This article is included in the Water and Global Change (WATCH) special collection.