Correction of Global Precipitation Products for Orographic Effects

Abstract

Underestimation of precipitation in topographically complex regions plagues most gauge-based gridded precipitation datasets. Gauge locations are usually in or near population centers, which tend to lie at low elevations relative to the surrounding terrain. For hydrologic modeling purposes, the resulting bias can result in serious underprediction of observed flows. A hydrologic water balance approach to develop a globally consistent correction for the underestimation of gridded precipitation in mountainous regions is described The adjustment is based on a combination of the catchment water balances and variations of the Budyko EIP versus PET/F curve. The method overlays streamflow measurements onto watershed boundaries and then performs watershed water balances to determine "true" precipitation. Rather than relying on a modeled runoff ratio, evaporation is estimated using the Budyko curves. The average correction ratios for each of 357 mountainous river basins worldwide are spatially distributed across the basins and are then interpolated to ungauged areas. Following application of adjustments for precipitation catch deficiencies, the correction ratios are used to scale monthly precipitation from an existing monthly global dataset (1979-99, 0.5° resolution). The correction for orographic effects resulted in a net increase in global terrestrial precipitation of 6.2% (20.2% in orographically influenced regions only) for the 1979-99 climatology. The approach developed here is applicable to any precipitation dataset in regions where good streamflow data exist. As a cautionary note, the correction factors are dataset dependent, and therefore the adjustments are strictly applicable only to the data from which they were derived. [PUBLICATION ABSTRACT]

Headnote

ABSTRACT

Underestimation of precipitation in topographically complex regions plagues most gauge-based gridded precipitation datasets. Gauge locations are usually in or near population centers, which tend to lie at low elevations relative to the surrounding terrain. For hydrologic modeling purposes, the resulting bias can result in serious underprediction of observed flows. A hydrologic water balance approach to develop a globally consistent correction for the underestimation of gridded precipitation in mountainous regions is described The adjustment is based on a combination of the catchment water balances and variations of the Budyko EIP versus PET/F curve. The method overlays streamflow measurements onto watershed boundaries and then performs watershed water balances to determine "true" precipitation. Rather than relying on a modeled runoff ratio, evaporation is estimated using the Budyko curves. The average correction ratios for each of 357 mountainous river basins worldwide are spatially distributed across the basins and are then interpolated to ungauged areas. Following application of adjustments for precipitation catch deficiencies, the correction ratios are used to scale monthly precipitation from an existing monthly global dataset (1979-99, 0.5° resolution). The correction for orographic effects resulted in a net increase in global terrestrial precipitation of 6.2% (20.2% in orographically influenced regions only) for the 1979-99 climatology. The approach developed here is applicable to any precipitation dataset in regions where good streamflow data exist. As a cautionary note, the correction factors are dataset dependent, and therefore the adjustments are strictly applicable only to the data from which they were derived.

(ProQuest Information and Learning: ... denotes formulae omitted.)

1. Introduction

Macroscale hydrological models are applicable to quantitative assessments of water and energy budgets on continental and global scales (e.g., Levis et al. 1996), to climate monitoring and analysis (e.g., Mote et al. 2005; Oltchev et al. 2002), and as a link between atmospheric and water resources models (e.g., Christensen et al. 2004). Precipitation is the main driver of the land surface hydrologic system and, therefore, hydrologic simulations are highly sensitive to misspecification of the precipitation forcings. Precipitation datasets that are constructed by the interpolation of point estimates to a coarse-resolution grid generally misrepresent (usually underestimate) precipitation in topographically complex regions due to an underrepresentation of gauge locations at high elevations. Therefore, interpolation algorithms that use simple weighting schemes like inverse distance are unable to capture the influence of orographie lifting on precipitation. Various publications (including Wallace and Hobbs 1977; Barros and Lettenmaier 1994; Roe 2005) discuss the influence of topography on precipitation and demonstrate that it depends on the mean flow field near the surface, that is, whether the slope is windward (upslope) or lee side (downslope) of the orographic barrier. In short, increased condensation and precipitation on the windward side and decreased precipitation on the leeward side result in a precipitation divide upwind of the orographic divide. The more arid downslope has been termed the "rain shadow." Roe (2005) reviews the mechanisms of orographie precipitation and the state of the art in observing and modeling of orographic precipitation.

A number of gridded (mostly monthly) precipitation time series are available on a global scale (e.g., Mitchell et al. 2005, manuscript submitted to J. Climate, hereafter M05; Chen et al. 2002; Willmott and Matsuura 2001; New et al. 2000; Huffman et al. 1997). These global precipitation datasets have been used for a number of purposes, including diagnoses of weather and climate models, evaluation of long-term changes in climate, forcing of hydrology and water resources models, and many other purposes. Most of these datasets have been derived from station data (and, in some cases, other data sources such as satellite algorithms and surface radars) that do not account for orographie effects. Using the Global Precipitation Climatology Project (GPCP) product (Huffman et al. 1997) as a driver for their simulations, Nijssen et al. (2001a) describe the difficulty in predicting streamflow for large river basins with considerable orography (in particular the Brahmaputra and Columbia River basins) and discuss the need for improved gridded precipitation estimates that account for orographic effects. One such product is based on the Parameter-Elevation Regressions on Independent Slopes Model (PRISM) (Daly et al. 1994, 2001, 2002). PRISM uses a digital elevation model (DEM) to partition the topography into facets of reasonably constant slope orientation. Precipitation for a particular grid cell is estimated by linearly regressing precipitation against DEM elevation for the stations on that cell's topographic facet. PRISM produces a monthly climatology for 1961-90 at a resolution of 1.25' for several countries, including western Canada (Alberta, British Columbia, Manitoba, Saskatchewan, and Yukon Territory), China, Mongolia, Puerto Rico, Taiwan, and the United States (Daly et al. 1994, 2001, 2002). Recent hydrological simulations for many large U.S. river basins have utilized PRISM to significantly improve streamflow predictions (e.g., Maurer et al. 2002). Although PRISM products have been developed for a few countries, they are not available globally, in part because of the model's requirement of a dense station network distributed over a range of elevations.

The objective of this study was to develop a method to account for orographie effects in 0.5° gridded precipitation estimates on a global scale. The method was designed to utilize existing gridded precipitation products, in particular the 1950-99 monthly time series of Willmott and Matsuura (2001). In previous work (Adam and Lettenmaier 2003), we developed adjustments for the systematic undercatch of precipitation by gauges. These adjustments, which accounted for the effects of wind-induced undercatch of liquid and solid precipitation and wetting losses on monthly means for a 0.5° latitude-longitude grid, were specific to the 1979-98 period. Application of these "catch ratios" to the Willmott and Matsuura (2001) data resulted in a net increase of 11.7% in global terrestrial mean annual precipitation (Adam and Lettenmaier 2003). The orographie adjustments developed in this study are specific to the same Willmott and Matsuura (2001) dataset and the slightly longer 1979-99 time period. Our final product is a 0.5° grid of correction ratios for all global land areas (excluding Antarctica and Greenland) that can be used to scale the mean annual precipitation. As a cautionary note, the correction ratios developed here are strictly applicable only to the data from which they were derived, and application of this method to other datasets could result in significantly different correction ratios. Note that hereafter the words "adjustment" and "correction" are used interchangeably, although strictly speaking one might argue that adjustment is more accurate because correction implies a priori knowledge of the final estimates.

Our approach was to distribute streamflow measurements over orographically influenced watersheds globally and to perform a water balance to estimate what the long-term precipitation should have been to produce the observed runoff. This approach is similar to one described by Fekete et al. (2000) in which runoff ratios were simulated on a global 0.5° grid using a simple water balance model. The authors then used these runoff ratios and their gridded estimates of runoff (which are a composite of simulated runoff and observed streamflow distributed onto the watershed) to calculate a new precipitation value. Given problems with this approach (e.g., an assumption that the simple water balance model is correctly simulating the spatial variation of runoff, which is not the case because the input precipitation data do not account for orographie effects), we opted not to apply it. Rather than relying on a modeled runoff ratio (and assuming it to be independent of precipitation depth over a basin), we used a variant of the Budyko E/P versus PET/P curve to compute the water balance. This relationship and the longterm water balance were solved simultaneously for mean annual basin average precipitation. This calculation, hereafter termed the "Budyko method," was the foundation of our approach, although it was only the first in a series of steps in creating a grid of 0.5° adjustment ratios.

2. Analytical methods

a. Overview of methodology

Because gridded precipitation estimates can be considered to be more accurate in regions of low topographic gradients (both because the effects of orography are less and because station locations tend to be more uniformly distributed where complex topography does not preclude access), we limited our analysis to areas with strong orographic influence. The preselection of a correction domain, which was based on topographic slopes derived from a DEM, is described in section 2b. Within this domain, we created a 5-min grid of correction bands that is a categorization of the degree of orographic influence for each grid cell (also described in section 2b). This categorization was based on relative elevation differences within a predefined radius of influence. Definition of the correction domain and creation of the correction band grid set the framework for the determination of the correction ratios, which was divided into three steps. In step 1 mean annual basin average precipitation was determined for 524 selected (357 mountainous and 167 low relief) river basins worldwide using the Budyko method (described in section 2c). Comparison of this "true" precipitation to the unadjusted precipitation yielded a basin-average correction ratio for each of these basins, which were then adjusted for bias using the low relief basins. Using the 5-min grid of correction bands, a spatial distribution of correction ratios was created for each of the gauged mountainous basins (step 2, described in section 2d). Finally, correction ratios in the gauged basins were interpolated over the entire correction domain using the grid of correction bands and knowledge of the dominant effective wind direction to separate upslope cells from downslope cells (step 3, described in section 2e). The analysis involved in both steps 2 and 3 was performed at a resolution of 5' (approximately 9 km × 9 km near the equator), a scale that arguably captures broad-scale topographic features (Daly et al. 1994). The final correction ratios were aggregated to 0.5° prior to application to the unadjusted precipitation grids.

b. Definition of correction domain

The DEM used to define the correction domain (and used throughout the study) is a global 30-arcsecond dataset derived from the Shuttle Radar Topography Mission (SRTM) (Farr and Kobrick 2000), a recent improvement to GTOPO30 (1996). This dataset is hereafter described as SRTM30. The 30-arcsecond data were aggregated to a 5' resolution prior to all subsequent analysis. We defined the initial correction domain as all 0.5° grid cells with slopes (aggregated from 5') greater than 6 m km^sup -1^. This threshold is the approximate slope above which the Willmott and Matsuura (2001) 1979-99 precipitation climatology differs by more than 10% from PRISM (Daly et al. 1994, 2001, 2002) over the contiguous United States. We determined this threshold by calculating the relative difference for each half-degree grid cell between the Willmott and Matsuura dataset and that of a half-degree aggregation of PRISM, binning the relative differences and calculating the mean of the terrain slopes for the cells in each bin. The final correction domain was created by running the initial correction domain through a spatial smoothing process using a radius of three 0.5° grid cells, that is, a grid cell was included in the correction domain if at least half of the cells within a threecell radius had a slope greater than threshold. The 0.5° slopes are shown in Fig. 1; the final correction domain is within the heavy black boundary lines.

All 5' grid cells within the correction domain were assigned to a correction band ranging from 2 (lowest elevations, least orographic influence) to 7 (highest elevations, most orographic influence). The assignment of correction bands was performed separately for each 0.5° cell in the correction domain. The bands were assigned by determining the maximum and minimum elevations of all 5' cells within a specified number of grid cells away from the overlying 0.5° cell and evenly dividing the elevations between minimum and maximum into the six correction bands. A radius of two 0.5° cells was used for the analysis, which corresponds to a radius of influence of approximately 125 km. The reasoning is that the degree of topographic influence is a function mostly of the relative differences in elevation and not as much of the absolute elevations. The 5' correction band assignments are shown in Fig. 2. All 5' cells outside the correction domain are assigned a value of one (not shown).

c. Determination of average correction ratios for gauged basins: Step 1

1) THE BUDYKO METHOD The continuity equation applied to a watershed is valid across all spatial and temporal scales:

... (1)

where P, E, and Q are the basin-average precipitation, evapotranspiration, and runoff, respectively; G is the net discharge of groundwater out of the aquifer underlying the basin; and dS is the net change in storage for a given time increment, dt. For longer time periods in which the net change in storage is negligible (e.g., reservoir and aquifer storage effects are not significant), Eq. (1) becomes

... (2)

where ... are long-term mean annual basinaverage precipitation, evapotranspiration, and runoff, respectively. Therefore, the precipitation climatology for a watershed can be determinedly distributing mean annual streamflow measurements ... back onto the watershed (i.e., dividing the streamflow volume by the area of the watershed).

Because evaporation is also an unknown, a second equation is needed. Building on the work of others (e.g., Schreiber 1904; Ol'dekop 1911), Budyko (1974) developed a relationship between the evaporative index, ..., and the aridity index, φ = ... (where ... is basin-average potential evapotranspiration). This relationship can be termed the Budyko ... versus ... curve:

.... (3)

The curve is semiempirical: the limits reflect physical constraints, but the curve was developed from a large number of observations. The Budyko curve and the physical constraints are shown in Fig. 3. The evaporative index is limited either by the energy available to evaporate moisture (in which ... approaches ...) or by the amount of moisture available for evaporation (in which ... approaches unity). After testing the relationship in 1200 regions, Budyko and Zubenok (1961) show that the mean discrepancy between the evaporative index as calculated from Eq. (3) and that derived by the water balance was about 10%. Budyko (1974) argues that this relationship can be satisfactorily applied to most mountainous basins, concluding that "the observed changes in total runoff and the runoff coeffitient are fully explained by the increase of precipitation, that is, by climatic factors," although he warns that this may not be universally true for the highest mountain basins or where there is a sharp change in the character of the underlying surface with increase of altitude. Furthermore, Budyko (1974) stresses that this relationship should only be applied to watersheds of "considerable size" where runoff does not vary appreciably under the influence of local conditions. Although Budyko did not define "considerable," his early testing of the relationship was only for basins exceeding 10 000 km^sup 2^ (Budyko 1951). Therefore, when performing our analysis, we attempted to select basins with drainage areas that exceeded 10 000 km^sup 2^.

Various improvements to the Budyko ... versus ... curve have been published in recent years. Milly (1994) incorporated several other variables into the estimation of the evaporative index in a rigorous mathematical framework. Zhang et al. (2001) incorporated the "plant-available water coefficient" (representing the relative difference in the way plants use soil moisture for transpiration) into the relationship. Sankarasubramanian and Vogel (2002) took into account the soil moisture storage capacity in their model and developed relationships for 458 watersheds throughout the United States with at least 10 yr of streamflow records. The Sankarasubramanian and Vogel (2002) ... versus ... curves for three values of γ, the soil moisture storage index, are shown in Fig. 3. We chose to use the Sankarasubramanian and Vogel (2002) relationship because of the practicality of applying it on a global scale and the unambiguity of input parameters [implementation of the Milly (1994) method is not practical on a global scale and the plant-available water coefficient of Zhang et al. (2001) was intended as a calibration parameter]. Furthermore, although the Sankarasubramanian and Vogel (2002) equation was developed only for U.S. river basins, basins from all of the contiguous U.S. climate types were incorporated, and these capture many of the climate types found worldwide. The evaporative index is a function of φ and γ as follows:

... (4)

where R = exp(-φ/γ) and γ = .... Sankarasubramanian and Vogel (2002) state that one possible parameterization for b is b^sub max^ = max(E^sub t^) + max(S^sub t^) in which max(S^sub t^) is the maximum soil moisture holding capacity of the basin, and max(E^sub t^), the maximum of the evapotranspiration state variable, can be found from the asymptotes of Budyko's framework:

... (5)

for which we used an annual time step, t. Assuming that ..., and max(S^sub t^) are known, the unknowns ... and ... were found by solving Eqs. (2) and (4), which was performed numerically. Finally, we calculated a basin-average correction ratio R^sub ave^ by dividing ... found using the Budyko Method by ... derived from the data of Adam and Lettenmaier (2003). Assuming that the gridded precipitation estimates are accurate outside of our correction domain (low relief regions), we determined if there was a consistent bias in the Budyko method by applying the method over low relief basins on each continent. We used this knowledge to adjust the bias in the Budyko basin-average correction ratios for the mountainous basins, which is discussed in section 2c(3).

2) PREPROCESSING OF DATA NEEDED FOR THE BUDYKO METHOD

Global Runoff Data Centre (GRDC: see online at http://grdc.bafg.de), River Discharge (RivDis) (Vörösmarty et al. 1998), and Hydro-Climatic Data Network (HCDN) (Slack and Landwehr 1992) datasets were utilized to select the river basins needed for the Budyko method. Figure 4 shows all of the streamflow gauging stations that have at least one year of observations between 1979 and 1999 and have a drainage area greater than 1000 km^sup 2^. A river basin was selected if its drainage area overlapped by at least 50% with the correction domain. We tried to choose river basins with drainage areas greater than 10 000 km^sup 2^ and with at least four years of record coincident with our period. However, in order to make the basin coverage greater, we relaxed these requirements in a few cases. Figure 5 shows the coverage of the 357 mountainous basins (Fig. 5a) and the 167 low relief basins (Fig. 5b) chosen for the Budyko method. There are large parts of the correction domain (gray areas in Fig. 5) that do not have coverage (e.g., parts of Africa, Asia, and South America), and the validity of our approach will be limited in those areas. Table 1 lists the number of basins selected for each continent and the percentage of those basins having less than four years of streamflow record. (Note that Borneo, New Guinea, and New Zealand are included in the "Australian" domain, whereas continental Australia does not contain any of the correction domain.) Figure 6 shows the distribution of basins (worldwide) according to streamflow record length. Table 2 lists the number of basins for each range of drainage area and shows that the majority of basins are greater than 10 000 km^sup 2^ (68% and 90% for mountainous and low relief, respectively). We used the basin delineations of the HYDRO1k digital dataset, which use the basin-coding scheme of Verdin and Verdin (1999). Mean annual basin runoff depth ... was determined for each of the basins by dividing the mean annual basin volume (for years coincident with 1979-99) by the drainage area of the basin.

The concept of PET was first introduced by Thornthwaite (1948) and was intended to be an index of the "drying power" of the climate. It is the rate at which evapotranspiration would occur from a surface completely covered with growing vegetation with an unlimited supply of soil moisture and without advection or heat storage effects (Dingman 2002). Because evapotranspiration is also dependent on vegetation characteristics, it is necessary to define a reference crop for PET. Allen et al. (1994) defined the reference crop as having a crop height of 0.12 m, a canopy resistance of 70 m s^sup -1^, and an albedo of 0.23, which is the reference crop assumed in this study. The FAO-56 Penman-Monteith equation (Allen et al. 1998) for estimating PET is widely used but has many data requirements that make application of the equation difficult in data-sparse areas. Other less data-intensive equations exist (with attendant assumptions), such as the temperature-based equation of Hargreaves and Samani (1982). Droogers and Allen (2002) make an argument that, given the inaccuracy of various meteorological observations in more remote parts of the globe, a modified version of the Hargreaves equation is best when limited climatological data are available or when inaccurate weather data collection is expected. This equation is similar to the Hargreaves equation but includes a precipitation term that acts as a surrogate for relative humidity:

PET = 0.0013S^sub 0^(T^sub avg^ + 17.0)(TD - 0.0123P)^sup 0.76^, (6)

where T^sub avg^ is mean temperature (°C) for a given day, TD is the difference between mean daily maximum and mean daily minimum temperature (°C) for a given month, P is the precipitation for a given month (mm), and S^sub 0^ is the water equivalent for extraterrestrial solar radiation (mm day^sup -1^) (Droogers and Allen 2002). The extraterrestrial solar radiation is given by

S^sub 0^ = 15.392d^sub r^(ω^sub s^ sinφ sinδ + cosφ cosδ sinω^sub s^), (7)

where d^sub r^ is the relative distance between the earth and the sun [Eq. (8)], ω^sub s^ is the sunset hour angle in radians [Eq. (9)], φ is the latitude in radians, δ is the solar declination in radians [Eq. (10)], and J is the Julian day number (Shuttleworth 1992):

... (8)

... (9)

... (10)

In polar regions, the sunset hour angle ω^sub s^ was set to zero when the sun never rises (i.e., when tanφ tanδ is less than -1) and to π when the sun never sets (i.e., when tanφ tanδ is greater than 1). PET was calculated for every day between 1979 and 1999 using the monthly temperature data of M05, the monthly precipitation data of Adam and Lettenmaier (2003), the daily temperature data of Sheffield et al. (2004) through 1995, and an extended version of the Nijssen et al. (2001b) data after 1995. Daily PET was aggregated to annual before application of the Budyko Method, and P... was determined for each of the 524 basins shown in Fig. 5 for the years coincident with the streamflow record. The 0.5° dataset of Dunne and Willmott (1996) was used for maximum soil moisture storage capacity, max(S^sub t^). A basin-average value was determined for each of the 524 basins shown in Fig. 5.

3) BIAS-ADJUSTMENT OF THE BUDYKO METHOD

By calculating the Budyko R^sub ave^ for the 167 low relief basins worldwide, we found that there was a consistent negative bias in the Budyko method, that is, the basin-average correction ratios were less than one for most of the basins. We found that there is a slight positive dependence between the aridity index φ and the bias, defined as the difference between R^sub ave^ and one (for the low relief basins). There is a priori expectation of this dependence because the need for streamflow diversion and storage and also the evaporative losses associated with irrigation and reservoir storage increase as φ increases, thereby causing a larger evaporative index with respect to φ than suggested by the Sankarasubramanian and Vogel (2002) formulation [Eq. (4)]. Therefore, using φ as a means to develop a weighting system, the bias for each mountainous basin on a particular continent was determined by taking a weighted average of the biases for low relief basins as follows:

... (11)

where bias, and W^sub i^ are the bias and weighting for the ith low relief basin, respectively, and n is the number of low relief basins on that continent. The φ ratios for both the low relief and the mountainous basins were sorted into 10 bins evenly divided between the minimum and maximum and assigned a bin number. The weight was calculated as

W^sub i^ = exp(10 - bindist^sub i^), (12)

where bindist^sub i^ is the number of units between the bin number for the mountainous basin and the bin number for the ith low relief basin. An exponential distribution was used in order to give very high weights to basins with similar φ ratios. The bias-adjusted basin-average correction ratio was calculated as

R^sub ave^ = R^sub ave,orig^ - bias, (13)

in which R^sub ave,orig^ is the value calculated using the Budyko method. For purposes of testing, this procedure was also applied to the low relief basins.

We calculated, for the low relief basins, the relative mean absolute error (RMAE), the square of the correlation coefficient (R^sup 2^), and the relative bias (RB) calculated as

... (14)

for basin-average precipitation between the gridded data of Adam and Lettenmaier (2003) (...) and the value predicted by the Budyko method (...) before and after bias adjustment (Table 3). (Note that there was no bias adjustment performed for the Australian correction domain, namely, the islands of Borneo, New Guinea, and New Zealand, owing to inadequate streamflow data.) After bias adjustment, the RMAE decreased, the R^sup 2^ increased, and the RB decreased to less than 4% for all continents, demonstrating that the bias adjustment was successful. Table 3 also shows, as an average for each continent, the resulting increases in precipitation due to bias adjustment for the mountainous basins. Figure 7 is a scatterplot of ... versus ... for eacn continent before and after bias adjustment (left and right panels, respectively). Linear regressions through the data for the low-relief basins approach the one-to-one line after bias adjustment. Note that the mountainous basins generally fall underneath the one-to-one line, demonstrating that the Budyko method predicts higher precipitation in the mountains.

d. Spatial distribution of correction ratios over gauged basins: Step 2

Spatial variability of the correction ratio across each of the gauged basins was constructed by developing a relationship between the correction ratio and the 5' correction bands (discussed in section 2b). PRISM (Daly et al. 1994, 2001, 2002) for "northwestern North America," defined hereafter as western Canada (Alberta, British Columbia, Manitoba, Saskatchewan, and the Yukon Territory), the contiguous United States, and Alaska, was used to determine the form of this relationship, and we therefore assumed that PRISM correctly captures the variability of precipitation with elevation. The PRISM climatology was aggregated from 1.25' to 5' for the analysis. The Budyko method was applied to determine basin-average precipitation ... and the basin-average correction ratio R^sub ave^ for 101 mountainous basins in northwestern North America. A set of scaling ratios (for correction bands 2 through 7) was determined for each basin such that, if applied to a 5' disaggregation of the 0.5° data of Adam and Lettenmaier (2003), a 5' precipitation grid result that has a basin average of P and the identical spatial variability as in PRISM. A plot of the correction band versus scaling ratio revealed that the data can be fitted with a quadratic relationship, for example,

r^sub band^ = Aband^sup 2^ + Bband + C, (15)

in which r^sub band^ is the correction ratio associated with each of the correction bands, and A, B, and C are parameters. A set of parameters was found for each of the 101 basins by regressing r^sub band^ against band and forcing r^sub band^ to equal one outside the correction domain (i.e., where band =1). The parameter A was found to have a slight dependence on R^sub ave^ and therefore was calculated as a function of R^sub ave^ via linear regression:

A = 0.015 - 0.022R^sub ave^. (16)

The other parameters were calculated by imposing two constraints on the relationship in Eq. (15). First, the correction ratio outside of the correction domain must be one (i.e., r^sub 1^ = 1), which gives

C = 1 - A - B. (17)

Second, R^sub ave^ must be conserved over the basin; that is,

... (18)

in which P^sub band^ is the sum of precipitation (of Adam and Lettenmaier 2003) for all 5' grid cells having the specified correction band. Substituting Eqs. (15) and (17) into Eq. (18) and solving for B yields

... (19)

Equations (15), (16), (17), and (19) result in a model in which the spatial variability of the correction ratio can be determined for all of the basins worldwide using ... and R^sub ave^ (both found using the Budyko method), the 5' correction band grid, and the precipitation to be adjusted (Adam and Lettenmaier 2003) disaggregated without interpolation to 5'. Implicit in this model is the assumption that the spatial variation of precipitation with correction band in this region is representative of this variation worldwide.

To examine how well the model recreated the original PRISM ratios, it was applied to each of the 101 northwestern North American basins and the RMAE was calculated between the original and modeled data points. Figure 8 is a histogram of the RMAE values, showing that 46% of the basins have RMAEs less than 10%, while 81% have RMAEs less than 20%. A further test was performed using PRISM data over Mongolia. Using five basins in Mongolia (ranging in size from 4600 to 445 000 km^sup 2^), the parameter A was found to vary with R^sub ave^ according to

A = -0.062 + 0.011R^sub ave^. (20)

The spatial variability model was applied to the Mongolian basins using both Eqs. (16) and (20), and the modeled r^sub band^ values were plotted against correction band for both cases (Fig. 9). The RMAE values for both cases are also shown in the figure. This test demonstrates that 1), although the scaling coefficients and offsets are considerably different in Eqs. (16) and (20) (and even opposite in sign), both parameterizations produce a reasonably good fit with low RMAEs, and 2) the quadratic form is adequate for Mongolian basins, indicating that the spatial variation model can be applied over continents other than North America, although ideally this model would be tested more rigorously for different regions if the means to do so existed.

e. Interpolation of correction ratios to ungauged basins: Step 3

Correction ratios were interpolated from grid cells in gauged basins to all grid cells in the correction domain using the correction band grid and a grid of slope types. (The correction ratios should be affected by the type of slope, e.g., the rain shadow occurs on the downslope, and the correction ratios for a specific band on a downslope will be different from those on an upslope.) All grids involved in the interpolation are at a 5' resolution. A linear distance weighting scheme was used in which the correction ratios were interpolated from grids with the same slope type and correction band. The weights, W^sub i^, were calculated as

... (21)

where D^sub i^ is the distance in kilometers between grid cells and R is the search radius in kilometers. The interpolated ratio was calculated as

... (22)

where n is the number of grid cells with matching slope type and correction band within the radius, R. A search radius of 500 km was used for the interpolation, but this radius was increased incrementally by 50 km if the minimum of 10 data points was not met. Finally, the 5' interpolated data were aggregated to 0.5°, and all 0.5° grid cells outside of the correction domain were set to one.

The 0.5° dominant effective wind direction was needed to determine slope type and was calculated using the National Centers for Environmental Research-National Center for Atmospheric Research (NCEP-NCAR) reanalysis (Kalnay et al. 1996) daily meridional and zonal wind speeds and precipitation between 1979 and 1999 (Fig. 10). For days that exceeded 0.5 mm in precipitation, the average daily wind direction was binned into one of the eight major wind directions, and the direction with the most occurrences for the 21-yr period was taken as the dominant effective wind direction for that grid cell. The direction of steepest slope (over a scale of approximately 50 km) was computed using the 5' DEM (aggregated SRTM30), and this 5' grid was overlaid with the 0.5° grid of dominant effective wind directions to determine slope type. Both wind direction and direction of steepest slope are in one of the eight major directions. If the wind direction (the direction the wind is traveling from) was the same as or was adjacent to the direction of steepest uphill, the slope was defined as downslope. Similarly, if the wind direction was the same as or was adjacent to the direction of steepest downhill, the slope was defined as upslope. The other possibilities were defined as cross slopes. The final slope type grid was created by running the calculated slope types through a spatial smoothing process using a radius of three 5' grid cells (Fig. 11).

3. Summary results

The 0.5° global correction ratios are shown in Fig. 12. These correction ratios were applied to the annual and monthly (1979-99) climatologies of the Adam and Lettenmaier (2003) data. The absolute increase in precipitation and the climatology after correction are shown in Fig. 13. Most of the corrections resulted in an increase in precipitation and, in some mountain ranges, the corrections exceeded a factor of 2, such as in the Andes (reaching a maximum of 3.7 in Peru), the Atlas Mountains (reaching a maximum of 2.5), the Zagros Mountains (reaching a maximum of 3.9), and the Himalayas (reaching a maximum of 3.5). In other areas, there was very little change in precipitation and even a decrease in precipitation (correction ratios less than one) such as in much of Europe, eastern Africa, isolated regions in Asia, and Borneo.

It is reasonable to expect low correction ratios on the lee side of major divides, but large isolated areas of low corrections suggest either that 1) the correction is unnecessary (in regions where precipitation does not change) or 2) there are deficiencies associated with the method. For example, we expect that orographic corrections in Europe are not as important as in other continents because the distribution of precipitation stations with elevation matches more closely the hypsometric curve. To evaluate this hypothesis, we obtained precipitation station information from the NCEP Climate Prediction Center (CPC) Summary of the Day data archived at NCAR. This dataset includes 15 190 stations globally and is representative of (although not identical to) the precipitation stations used to create gridded datasets; for example, Willmott and Matsuura (2001) utilized 20 782 stations and New et al. (1999) utilized 19 295 stations. Figure 14 shows the difference between the percent of stations and the percent of area for each elevation band globally using the Summary of the Day data; the figure confirms that the differences are lowest in Europe.

Outside of Europe, low correction ratios are most likely due to limitations in the methodology, such as the use of unnaturalized streamflow data; for example, river basins with significant diversions result in Budyko P values that are too low and therefore correction ratios that are too low. As discussed in section 5, much of Africa and eastern Asia are data-poor areas, and the Budyko method may not produce realistic results, even if naturalized streamflow data were used. This absence of bias adjustments over the Australian correction domain also partly explains the low correction ratios over Borneo.

The mean increase in precipitation due to the application of the correction ratios was computed globally and for each of the continents (Table 4). Eurasia and North America have the highest increases partly because a greater portion of these continents is affected by orographic processes. Table 4 also lists the increases within the correction domain only, for each continent. Increases were calculated for each elevation band (in increments of 500 m) for all continents (Fig. 15) using the 5' aggregation of SRTM30. In general, the increase in precipitation due to orographic correction increases with elevation (with the exception of Africa), although mean Eurasian corrections reach a maximum between 2 and 3.5 km.

4. Analysis and discussion

a. Effects of interpolation: Conservation of magnitude and spatial variability

Although P for each basin found using the Budyko method (step 1) was forced to be conserved for each basin when the spatial variability was created (step 2), it was not necessarily conserved during the interpolation process (step 3). To examine the extent of the effect of the interpolation process, plots of basin-average precipitation before and after correction versus the Budyko P were plotted for each continent (Fig. 16). Linear regressions were fitted through both the uncorrected and the corrected data, and R^sup 2^ values were calculated between the two sets of data and the Budyko P values (all shown in Fig. 16). In every case the R^sup 2^ values (initially between 6% and 87%) improved (to between 79% and 98%). In all cases the slope of the linear regressions improved; that is, they approached the one-to-one line.

For the 101 North American basins used to develop the spatial variation of the correction ratio within the gauged basins [Eqs. (15), (16), (17), and (19) in section 2d], we examined the degree to which the spatial variation of the correction ratios was conserved through each step. Figure 17a shows the correction ratios needed to create the Budyko P with the spatial variability of PRISM over northwestern North America. These are the correction ratios used to develop the PRISM-derived model for the spatial variation of the correction ratios across the gauged basins. Figure 17b shows the correction ratios for the same basins using the PRISM-derived model (step 2). Figure 17c shows the correction ratios after the interpolation step (step 3) but prior to the 0.5° aggregation. After interpolation, it is apparent that the variability of the correction ratios was reduced. The coefficients of variation (C^sub v^) of the correction ratios are 0.24, 0.23, and 0.15 for Figs. 17a, 17b, and 17c, respectively. The RMAE of the correction ratios is 9.0% between Figs. 17a and 17b, 9.7% between Figs. 17b and 17c, and 13.7% between Figs. 17a and 17c. The RB is 0.9% between Figs. 17a and 17b, 0.4% between Figs. 17b and 17c, and 0.5% between Figs. 17a and 17c. This suggests that, although the variability became smoother after interpolation, the aggregate effect of this smoothing was small for these 101 basins.

b. Comparisons with PRISM climatology

A 5' aggregation of the PRISM climatology over northwestern North America was used to derive the spatial variability of the correction ratios, but the magnitude of the corrections was determined using the Budyko method, which is independent of PRISM. Therefore, the PRISM climatology is an independent estimate of precipitation magnitude and provides a comparison for our corrected precipitation data. There are two challenges to this comparison. First, the climatology period is different; the PRISM climatology period is 1961-90, while our period is 1979-99. We made no attempt to account for this difference in the following comparisons. Second, in addition to orographic corrections, our precipitation data have been adjusted for gauge undercatch biases (Adam and Lettenmaier 2003), whereas the PRISM precipitation data have not. Therefore, rather than comparing precipitation depths, we compared the increase in precipitation due to our orographic corrections to the increase inferred by PRISM by using the Willmott and Matsuura (2001) 0.5° data as the baseline precipitation (the PRISM climatology was first aggregated to 0.5°). Our orographic corrections resulted in a net increase of 16.1% over northwestern North America and a net increase of 41.6% within the correction domain only. The PRISM-inferred corrections resulted in a net increase of 13.6% over northwestern North America, a net increase of 35.5% over the correction domain only and a net decrease of 0.3% outside of the correction domain. This yields a relative bias of 17% for the entire region and 18% over the correction domain.

Figure 18a shows the spatial distribution of the increases inferred by the PRISM corrections over the common region, and Fig. 18b shows the spatial distribution of the increases due to our orographic corrections. The correction domain is shown within the purple boundary lines. Although the broad features of the spatial variability are similar (which is expected because PRISM was used to aid in the development of this variation), the variation in Fig. 18b is smoother and the extremes are greatly reduced. (Within the correction domain only, the C^sub v^s of the increases are 1.41 and 0.42 for Figs. 18a and 18b, respectively, resulting in a 70% decrease in variability.) The smoothing was not due to parameterization of the variability (i.e., there is little decrease in variability between Figs. 17a and 17b), but rather to two aspects of the interpolation procedure. First, to smooth discontinuities at basin boundaries and to maximize the information available during interpolation, we used a search radius of 500 km (although the values at these distances had lower weights). The resuiting decrease in variability (about 35%) can be seen through comparison of Figs. 17b and 17c. Second, the interpolation procedure was accomplished using a weighted average rather than a functional relationship (such as regression) and, therefore, correction ratios for ungauged regions that would have extremely high or low corrections were sampled from data that may not represent extremes.

As a further comparison, the increases corresponding to Figs. 18a and 18b were determined for each elevation band in increments of 500 m (Fig. 19). Both datasets result in precipitation corrections that increase with elevation, although the shape of this increase is slightly different. As a final comparison, the precipitation increases corresponding to the data in Figs. 18a and 18b were separated by month (Fig. 20). Our orographic corrections are on a mean annual basis; therefore the variation of the increases by month is only dependent on the monthly precipitation prior to adjustment. On the other hand, the PRISM method is performed individually for each month, which results in a more realistic monthly variation; that is, the corrections are smaller for summer precipitation over northwestern North America. This limitation of our methodology is discussed further in section 5.

5. Limitations

a. Streamflow data

In some areas, there were few or no steamflow gauging stations that met our criteria for record lengths. Therefore, correction ratios in these areas were interpolated from stations far away, possibly from locations in regions that are not climatologically similar. Areas most affected by this limitation are central and northern China, southern Mongolia, the Middle East, the central Andes, northeastern Russia, and the western coast of Africa. Furthermore, in order for the Budyko method to be accurate, river basins need to be large enough for river runoff not to be appreciably influenced by local conditions of a nonclimatic nature (Budyko 1974), which we took to be larger than approximately 10 000 km^sup 2^. On the other hand, larger river basins may experience more diversions, which would result in correction ratios that are too low. We have made no attempt to naturalize the streamflow or to account for streamflow diversions, as the data that would be required are not generally available. As discussed in section 3, this could be a possible explanation for why correction ratios are uniformly low in several regions, for example, most of Africa and parts of Asia (although the bias adjustment procedure corrects for much of the low bias). Also, as mentioned previously, 16% of the stations globally have fewer than four years of record coincident with our period of interest. Basin storage effects (including reservoirs) may not be negligible for this short period, and this could cause misestimation of the correction ratios.

b. PET

Using the 33 mountainous basins in Africa, we performed a sensitivity study to see how the basin-average correction ratios (R^sub ave^) vary with changes in basin-average PET, runoff, and soil moisture storage capacity. The resulting curves are shown in Fig. 21. For the African basins, PET is the most sensitive of the three variables. [Because PET values are relatively high for the African basins, the formulation for the parameter b causes the soil moisture capacity to have little effect (see section 2c). In places with less PET, this would not be as much the case.] The reliance of this approach on PET and the high degree of dependence of our results on PET can be considered a limitation. If there is a net negative bias in the estimation of PET via the modified Hargreaves method (Droogers and Alien 2002) over Africa, the correction ratios would be biased downward, which could be a partial contributor to the poor results there.

c. Seasonal and interannual variability

Our methodology is a correction to the mean annual precipitation climatology. The Budyko method assumes negligible storage effects (which is approximately true on a long-term basis for which there are no significant streamflow diversions and little reservoir evaporation). Therefore, this method cannot be applied to estimating precipitation on a seasonal or interannual basis. The seasonal and interannual variability of the corrected precipitation data based on our method will therefore be the same as that of the unadjusted data. As mentioned in section 4, the effects of orography on precipitation usually are seasonal, for example, PRISM-inferred corrections over northwestern North America are smallest in May (24%) and greatest in October (53%) (Fig. 20). To investigate the error due to applying an annual adjustment, we derived a PRISM-inferred annual correction ratio by dividing the PRISM annual climatology by the Willmott and Matsuura (2001) annual climatology for each 0.5° grid cell in northwestern North America. Multiplying the Willmott and Matsuura (2001) monthly climatology by these ratios results in a dataset that has the same annual climatology of PRISM but the monthly variability of Willmott and Matsuura (2001). Finally, we averaged both the PRISM monthly climatology and the annually adjusted Willmott and Matsuura (2001) monthly climatology over the entire domain and over the correction domain only and calculated the relative differences with respect to PRISM (Table 5). The differences were all less than 10%, and so an annual (versus monthly) orographic correction is not a major limitation, at least at the continental scale.

6. Recommendations

Although our goal was to develop a spatially consistent method to correct for orographic effects on a global scale, we acknowledge that many potential applications of the resulting data are regional. Our approach is most appropriate for regions with high-quality long-term streamflow data for basins with few diversions and long-term storage effects, and over regions that are not heavily developed (i.e., where changes in the hydrologic cycle are mainly governed by climatic and not anthropogenic effects). Furthermore, the coverage of river basins should be adequate to fully represent the climate and topography found throughout the region. Therefore, a satisfactory degree of coverage will depend on the region of interest, and we leave this to the insight of the user. Table 6 gives the percentage of correction domain covered by basins for 22 regions of the globe and may be helpful to users who wish to apply these data in any one of these regions.

7. Summary

Our work is intended to satisfy a need for gridded precipitation data that account for orographic effects in a globally consistent framework (Nijssen et al. 2001a). Combined with the work of Adam and Lettenmaier (2003), the final product is a (1979-99) precipitation climatology that is adjusted for gauge catch biases on a monthly basis and orographic effects on an annual basis. Both adjustments are designed to utilize an existing 0.5° precipitation product (e.g., M05; Chen et al. 2002; Willmott and Matsuura 2001; New et al. 2000), but the grid of correction ratios for the orographic adjustment is specific to the dataset to be adjusted, that is, the Adam and Lettenmaier (2003) data, a derivative of Willmott and Matsuura (2001). Combination of the gauge catch deficiency and orographic adjustments resulted in a net increase of 17.9% (11.7% and 6.2%, respectively) of global terrestrial (excluding Antarctica) mean annual precipitation.

Acknowledgments. This research was supported by NASA Grant NAG5-9416 to the University of Washington.

References

REFERENCES

Adam, J. C., and D. P. Lettenmaier, 2003: Adjustment of global gridded precipitation for systematic bias. J. Geophys. Res., 108, 4257, doi:10.1029/2002JD002499.

Allen, R. G., M. Smith, A. Perrier, and L. S. Pereira, 1994: An update for the definition of reference evapotranspiration. ICID Bull, 43 (2), 1-34.

_____, L. S. Pereira, D. Raes, and M. Smith, 1998: Crop evapotranspiration: Guidelines for computing crop requirements. Irrigation and Drainage Paper 56, FAO, Rome, Italy, 300 pp.

Barros, A. P., and D. P. Lettenmaier, 1994: Dynamic modelling of orographically induced precipitation. Rev. Geophys., 32, 265-284.

Budyko, M. I., 1951: On climatic factors of runoff (in Russian). Probl. Fiz. Geogr., 16, 41-48.

_____, 1974: Climate and Life. Academic Press, 508 pp.

_____, and L. I. Zubenok, 1961: The determination of evaporation from the land surface. Izv. Akad. Nauk SSSR Ser. Geogr., 6, 3-17.

Chen, M. P., P. Xie, J. E. Janowiak, and P. A. Arkin, 2002: Global land precipitation: A 50-yr monthly analysis based on gauge observations. J. Hydrometeor., 3, 249-266.

Christensen, N. S., A. W. Wood, N. Voisin, D. P. Lettenmaier, and R. N. Palmer, 2004: Effects of climate change on the hydrology and water resources of the Colorado river basin. Climate Change, 62 (1-3), 337-363.

Daly, C., R. P. Neilson, and D. L. Phillips, 1994: A statistical-topographic model for mapping climatological precipitation over mountainous terrain. J. Appl. Meteor., 33, 140-158.

_____, G. H. Taylor, W. P. Gibson, T. W. Parzybok, G. L. Johnson, and P. Pasteris, 2001: High-quality spatial climate data sets for the United States and beyond. Trans. ASAE, 43, 1957-1962.

_____, W. P. Gibson, G. H. Taylor, G. L. Johnson, and P. Pasteris, 2002: A knowledge-based approach to the statistical mapping of climate. Int. J. Climatol, 16, 841-859.

Dingman, S. L., 2002: Physical Hydrology. 2d ed. Prentice Hall, 646 pp.

Droogers, P., and R. G. Allen, 2002: Estimating reference evapotranspiration under inaccurate data conditions. Irrig. Drain. Syst, 16, 33-45.

Dunne, K. A., and C. J. Willmott, 1996: Global distribution of plant-extractable water capacity of soil. Int. J. Climatol., 16, 33-45.

Farr, T. G., and M. Kobrick, 2000: Shuttle radar topography mission produces a wealth of data. Eos, Trans. Amer. Geophys. union, 81, 583-585.

Fekete, B. M., C. J. Vörösmarty, and W. Grabs, 2000: Global, composite runoff fields based on observed river discharge and simulated water balances. Documentation for UNH-GRDC Composite Runoff Fields, v.1.0, Global Runoff Data Centre, Koblenz, Germany, 114 pp.

GTOPO30, cited 1996: Global 30 arcsecond elevation data set. [Available online at http://edcwww.cr.usgs.gov/landdaac/ gtopo30/gtopo30.html.]

Hargreaves, G. H., and Z. A. Samani, 1982: Estimating potential evapotranspiration. J. Irrig. Drain. Eng., 108, 225-230.

Huffman, G. J., R. F. Adler, P. Arkin, A. Chang, R. Ferraro, A. Gruber, J. Janowiak, and A. McNab, 1997: The Global Precipitation Climatology Project (GPCP) combined precipitation dataset. Bull Amer. Meteor. Soc., 78, 5-20.

Kalnay, E., and Coauthors, 1996: The NCEP/NCAR 40-Year Reanalysis Project. Bull. Amer. Meteor. Soc., 77, 437-472.

Levis, S., M. T. Coe, and J. A. Foley, 1996: Hydrologic budget of a land surface model: A global application. J. Geophys. Res., 101 (D12), 16 921-16 930.

Maurer, E. P., A. W. Wood, J. C. Adam, D. P. Lettenmaier, and B. Nijssen, 2002: A long-term hydrologically based dataset of land surface fluxes and states for the conterminous United States. J. Climate, 15, 3237-3251.

Milly, P. C., 1994: Climate, soil water storage, and the average annual water balance. Water Resour. Res., 30, 2413-2156.

Mote, P. W., A. F. Hamlet, M. P. Clark, and D. P. Lettenmaier, 2005: Declining mountain snowpack in western North America. Bull. Amer. Meteor. Soc., 86, 39-49.

New, M. G., M. Hulme, and P. D. Jones, 1999: Representing twentieth-century space-time climate variability. Part I: Development of a 1961-90 mean monthly terrestrial climatology. J. Climate, 12, 829-856.

_____, _____, and _____, 2000: Representing twentieth-century space-time climate variability. Part II: Development of 1961-90 monthly grids of terrestrial surface climate. J. Climate, 13, 2217-2238.

Nijssen, B., G. O'Donnell, D. P. Lettenmaier, D. Lohmann, and E. F. Wood, 2001a: Predicting the discharge of global rivers. J. Climate, 14, 3307-3323.

_____, R. Schnur, and D. P. Lettenmaier, 2001b: Global retrospective estimation of soil moisture using the VIC land surface model, 1980-1993. J. Climate, 14, 1790-1808.

Ol'dekop, E. M., 1911: On evaporation from the surface of river basins (in Russian). Transactions on Meteorological Observations, University of Tartu, 4.

Oltchev, A., and Coauthors, 2002: The response of the water fluxes of the boreal forest region at the Volga's source area to climatic and land-use changes. Phys. Chem. Earth, 27, 675-690.

Roe, G. H., 2005: Orographic precipitation. Annu. Rev. Planet. Sci., 33, 1-31.

Sankarasubramanian, A., and R. M. Vogel, 2002: Annual hydroclimatology of the United States. Water Resour. Res., 38, 1083-1094.

Schreiber, P., 1904: Ueber die Beziehungen zwischen dem Niederschlag und der Wasseruhrung der Flusse in Mitteleuropa. Meteor. Z., 21, 441-452.

Sheffield, J., A. D. Ziegler, E. F. Wood, and Y. Chen, 2004: Correction of the high-latitude rain day anomaly in the NCEP-NCAR reanalysis for land surface hydrological modeling. J. Climate, 17, 3814-3828.

Shuttle-worth, W. J., 1992: Evaporation. Handbook of Hydrology, D. R. Maidment, Ed., McGraw-Hill, 4.1-4.53.

Slack, J. R., and J. M. Landwehr, 1992: Hydro-climatic data network (HCDN): A U.S. geological survey streamflow dataset for the United States for the study of climate variations, 1874-1988. U.S. Geological Survey Open-File Rep. 92-129, 193 pp.

Thornthwaite, C. W., 1948: An approach toward a rational classification of climate. Geogr. Rev., 38, 55-94.

Verdin, K. L., and J. P. Verdin, 1999: A topological system for delineation and codication of the Earth's river basins. J. Hydrol., 218 (1-2), 225-230.

Vörösmarty, C. J., B. Fekete, and B. A. Tucker, cited 1998: River discharge database, version 1.1 (RivDIS v1.0 supplement). [Available online at http://www.rivdis.sr.unh.edu/.]

Wallace, J. M., and P. V. Hobbs, 1977: Atmospheric Science: An Introductory Survey. Academic Press, 467 pp.

Willmott, C. J., and K. Matsuura, cited 2001: Terrestrial air temperature and precipitation: monthly and annual time series (1950-1999) (version 1.02). Center for Climate Research, University of Delaware, Newark, DE. [Available online at http://climate.geog.udel.edu/~climate/html_pages/archive. html.]

Zhang, L., W. R. Dawes, and G. R. Walker, 2001: Response of mean annual evapotranspiration to vegetation changes at catchment scale. Water Resour. Res., 37, 701-708.