Evapotranspiration and runoff from large land areas: Land surface hydrology for atmospheric general circulation models

A land surface hydrology parameterization for use in atmospheric GCMs is presented. The parameterization incorporates subgrid scale variability in topography, soils, soil moisture and precipitation. The framework of the model is the statistical distribution of a topography-soils index, which controls the local water balance fluxes, and is therefore taken to represent the large land area. Spatially variable water balance fluxes are integrated with respect to the topography-soils index to yield our large scale parameterizations: water balance calculations are performed for a number of intervals of the topography-soils distribution, and interval responses are weighted by the probability of occurrence of the interval. Grid square averaged land surface fluxes result. The model functions independently as a macroscale water balance model. Runoff ratio and evapotranspiration efficiency parameterizations are derived and are shown to depend on the spatial variability of the above mentioned properties and processes, as well at the dynamics of land surface-atmosphere interactions.


Introduction
One problem that climate modelers and hydrologists have in common is modeling the hydrologic cycle at large scales. Hydrologists have traditionally been interested in short and long term predictions of the water balance fluxes at the catchment scale (e.g. for flood forecasting and water supply predictions). However, there is no consensus on how to extrapolate well known point relationships to the catchment scale and beyond. Recently, hydrologists have also realized that processes operating on scales greater than the size of a watershed can be responsible for observed watershed response. For example, precipitation systems typically operate on scales larger than the catchment, and determining the origin of that water falling as precipitation may be a global scale problem (Koster, 1988). In short, hydrologists need to improve their understanding of hydrologic response at the large scale.
Climate modelers are particularly interested in large scale parameterizations of land surface hydrology for their numerical simulations of climate using atmospheric general circulation models (GCMs). But climate modelers face a scale problem of their own. Recognizing the sensitivity of GCM climate to land surface boundary conditions (e.g. Shukla and Mintz, 1982), climate modelers are seeking a more detailed grid scale parameterization than the simplified onedimensional approach first introduced into GCMs by Manabe (1969). Although the parameterization of land surface hydrological processes at the large scale is an area of active research (e.g. Dickinson, 1984;Sellers et al., 1986;Abramopolous et al., 1988;Entekhabi and Eagleson, 1989), the problem is still largely unresolved.
The hydrologic fluxes produced by large land areas are a result of the complex interaction of climate, soils, topography, geology, and vegetation. Each of these components exhibits a high degree of spatial variability at the large scale, and individual correlation lengths are not necessarily the same. For example, soil properties may vary widely on the hillslope scale, while vegetation and climate vary on more regional scales, Remote sensing can certainly help to quantify some of this variability, but the problem that remains is how to parameterize the dynamic hydrologic response of large land areas which exhibit considerable subgrid scale spatial variability in the major process controls.
The modeling problem is compounded when one considers that hydrologic phenomena at the land surface (i.e. infiltration, runoff, subsurface flow, and evapotranspiration) operate on different space-time scales. Infiltration and runoff, for example, are storm phenomena, which operate on the same time scales as precipitation events and on localized spatial scales. Subsurface flow, however, is active during storms and interstorm periods (as subsurface storm flow and base flow, respectively), and may act over the entire area of interest.
The scientific problem at hand is then to formulate mathematical descriptions for hydrologic response at the large scale, reconciling the various space-time scales of hydrologic phenomena with the tremendous spatial variability in land surface characteristics. Unfortunately, however, much of our knowledge of these processes is derived from laboratory experiments or field work at the point scale. Consequently, we must develop mathematical expressions for the large scale response rooted in the physics known at smaller scales, while incorporating important spatial variability in land surface hydrologic processes and properties.
A major control on the land surface hydrologic fluxes at the GCM grid square scale is subgrid scale variation in surface soil moisture. The dependence of local infiltration and evapotranspiration capacities on local surface soil moisture is well known. Consequently, as the framework of our parameterization, we have chosen to model the subgrid scale space-time dynamics of surface soil moisture. Given a local value of surface soil moisture, local infiltration, evapotranspiration, and surface runoff fluxes can be computed. Given the space-time distribution of surface soil moisture, local fluxes can then be aggregated to yield the areal averaged water balanced fluxes.
It is our belief that subgrid scale variations in topography and soil dominate the process of spatial redistribution of soil moisture over large land areas. In this paper we propose the use of hydrogeomorphic relationships, which incorporate these variations by means of a local topography-soils index (Beven, 1986), to model the space-time dynamics of soil moisture redistribution. Local land surface hydrologic fluxes can be related to surface soil moisture, and through the control of topography and soils on the local water balance, to local topography-soils indices. Therefore, knowledge of the distribution of the topography-soils index provides a consistent framework for averaging local land surfaces fluxes up to the macroscale.
Based on these concepts, we present a land surface hydrology parameterization for the areal average water balance fluxes, which includes a subgrid scale variability in precipitation, topography, soils, soil moisture, infiltration, evapotranspiration, and surface runoff. We have attempted to keep the model as simple as possible, while remaining true to the process physics, for two reasons: (i) to minimize computer run times, so that when combined with a more complicated GCM, no significant additional computational burdens arise, and (2) due to the lack of large scale data (e.g. soils, vegetation, evaporation) and enough variability to develop and verify more complex models. The model can be used by hydrologists to study the water balance over large land areas, and by climate modelers, to compute boundary conditions at the land surface for simulations of the general circulation of the atmosphere. Please note that throughout this paper, the term land area is used interchangeably with GCM land surface grid square, as the model is meant to meet the needs of both hydrologists and climatologists.
In the sections that follow, the dynamic soil moisture model is presented after a review of runoff generation mechanisms. Next, the equations for the land surface fluxes of infiltration, evaporation and transpiration are presented, with a description of how these flux rates are coupled to both surface soil moisture and atmospheric forcing. A discussion of runoff computations within the model is followed by a description of the spatial variability in rainfall. Finally, the macroscale equations for the land surface hydrologic fluxes are described. Horton (1933) first proposed that overland flow is generated over an entire catchment, whenever the infiltration capacity of its soil is exceeded by the precipitation rate. Today we realize that the concept of an entire watershed producing Hortonian overland flow cannot fully explain the generation of storm runoff. Numerous field studies have shown that at any time during a storm, only partial areas of the catchment are contributing runoff. Additionally, the spatial variability in soil properties, antecedent surface soil moisture, topography, and rainfall will result in a range of runoff generation mechanisms operating to produce the storm hydrograph. The common feature of these mechanisms is that the partial areas which contribute runoff are dynamic; they can expand and contract during storm and interstorm periods, and they can vary seasonally as well. The difference between these partial area mechanisms is the manner in which runoff is generated, and the pathway the runoff takes to the stream.

Runoff Generation Mechanisms
Hortonian overland flow, or infiltration excess runoff, occurs over those areas of a catchment where the local infiltration capacity is exceeded by the local precipitation rate. This is not a common occurrence in most humid, temperate climates. However, it may occur more frequently in arid and semi-arid climates, or where land use practices have severely altered the structure of the soil.
Rain falling directly on saturated soil areas or stream channels produces saturation excess runoff (Dunne and Black, 1970). These saturated areas are commonly located adjacent to streams, where the water table is high and soil moisture storage capacity is low. Saturation excess contributing areas grow during storms as infiltrating rainfall raises the shallow water table to the surface. The downslope redistribution of infiltrated rainfall toward the stream network also contributes to the growth of saturated areas along channels. The saturation excess runoff mechanism is most active in humid regions, with thin soils, gentle terrain, concave lower slopes, and wide valley bottoms (Dunne, 1978).
Subsurface storm flow occurs when infiltrated rainfall travels rapidly downslope, either through a network of interconnected large pores (macropores), or by flow in the saturated zone, as the water table rises during storms. Subsurface storm flow is most common in humid regions with steep, straight or convex slopes, narrow valley bottoms, and where permeable soils overlie relatively impermeable soils or bedrock (Dunne, 1978). Even when subsurface flow velocities are too small to contribute significantly to the storm by hydrograph, it may still dominate the overall response in terms of volume, as the delayed subsurface response provides the hydrograph tail (Knisel, 1973).
These runoff generation mechanisms are commonly in operation in many areas of the world. The dominant mechanism will be a function of climate and local land surface properties. However, any model for the generation of runoff from large land areas should incorporate these concepts.

A Physically Based Model for Space-Time Soil Moisture Dynamics
To properly model the land surface hydrologic fluxes, the spatial distribution of surface soil moistue becomes crucial. Near surface soil moisture influences infiltration and exfiltration capacities and thus the quantity of infiltration excess runoff and evapotranspiration. The amount of soil moisture storage available at any location determines the amount of saturation excess runoff produced. Subsurface flow is also dependent upon the moisture status of the area, contributing more as average soil moisture levels increase and less as they decrease. Beven and Kirkby (1979) present a simple, physically based model for the spatial distribution of soil moisture, based on topographic controls. The model was first proposed for small to medium sized catchments in humid temperature areas, but subsequent revisions (e.g. Sivapalan et al., 1987) allow for a much broader range of usage. The treatment of evapotranspiration in these papers, if present, has been weak, ignoring any spatial variations (Beven and Kirkby, 1979;Beven et al., 1984). In this paper we extend the concepts of Beven (1986) and Sivapalan et al. (1987) to large land areas, and we have included space-time variations in precipitation and evapotranspiration.

SOIL MOISTURE
To simplify the modeling process the following assumptions regarding land surface hydrologic prgcesses are made: (a) that the hydrologic response of large land areas proceeds as a series of quasi-steady states; (b) that the subsurface flow rate, qi, at location i can be related to a local storage deficit, si, by where Ti exp( -si/m) is conceptually equivalent to the local transmissivity, tan fl is the local slope, and pr is a parameter related to the decline of saturated hydraulic conductivity with depth (storage deficit in this paper refers to available porosity ' / beyond the field capacity of the soil -positive deficits imply unsaturated areas while negative deficits imply saturated areas); (c) that the saturated hydraulic conductivity declines exponentially with depth, so that where K, is the saturated hydraulic conductivity, z is the depth into the soil profile, K0 is the saturated hydraulic conductivity at the soil surface, and f is a scaling parameter; and (d) that at any location i, at steady state, the downslope saturated flow is given by where a is the unslope contributing area which drains through the unit contour width at i, and R is defined as the steady rate of recharge to the water table.
Integrating Equation (2) from a depth Z to the water table at depth zi, the transmissivity of the saturated soil profile, T(zi), is obtained. For large for Z where To = Ko/f Assuming a simple relationship between storage deficit and water table depth, e.g.
where A0 is the difference between saturation and residual moisture contents, and 6 is a constant, then (1) is equivalent to the downslope saturated flow beneath the water table when T~ = To, and m = 6AO/f, i.e. qi = To exp(-fzi)tan ft.
Note that these assumptions are not overly restrictive. Assumption (a) is a pragmatic decision given the large areas under consideration and the simplicity desired for compatibility with GCMs. Assumption (b) is nothing more than the kinematic wave assumption if we assume that the soil surface is roughly parallel to the bedrock and if the assumptions of the last paragraph are considered. The applicability of assumption (c) has been demonstrated by Beven (1982) for a wide variety of soils. Assuming a steady R in (d) is not unreasonable given no way to measure spatial variations in recharge rates on upslope contributing areas.
Equating (1)  In Te = ~ In To dA, and In(aTe~To tan fl) is the local value of the topography-soils index.
Equation (8a) states that the deviation of the local storage deficit from its areal average results from the deviation of the local value of In (a/tan fl) from its areal average and the deviation of the local transmissivity coefficient from its areal average. Restated, given any average moisture level g, the local storage deficit at steady state is determined by a topographic effect and a soils effect. Equation (8b) states that at any point on the land surface, for a given #, one need only the local value of the topography-soils index to determine the storage deficit.
To infer surface soil moisture from storage deficit the following assumptions are made. When the storage deficit is zero or negative, the soil is saturated and the surface soil moisture is equivalent to the saturation moisture content of the soil. When the storage deficit is positive, the location is unsaturated, and the soil profile is assumed to have reached gravity drainage. Local surface soil moisture content can then be determined from the Brooks and Corey (1964) Substituting (6) into (7), while solving (6) for R and using in (7), after some manipulation, an expression for storage deficit at any loction i is given as: and 0 is the moisture content, ~ is the matric head equal to zi -z, 0r is the residual moisture content, 0s is the saturation moisture content, ~b c is the height of the capiIlary fringe, and B is the pore size distribution index. At the soil surface, where z --0 and ~ = zi, the surface soil moisture, Oi, is The relationship between zi and s/, i.e. zi = si/~AO, is obtained by integrating the soil moisture characteristic from the top of the capillary fringe to the soil surface (Sivapalan et al., 1987). Therefore, for any unsaturated location, the water table depth and the surface soil moisture can be computed.
The spatial distribution of storage deficits and surface soil moisture can now be computed given g, the topography-soils index, and some soil parameters. The topography index can be calculated using digital elevation models (DEMs) or topographic maps. A number of algorithms are currently available to extract information such as upslope contributing area and slope from DEMs. Where available, data on hydraulic conductivity may be used to determine the topographysoils index. When these data are not available, the mean and standard deviation of hydraulic conductivity will suffice (discussed later). Also discussed in a later section is the estimation of g and f or m. For simplicity, the soil parameters Or, Os, d/c, and B are taken as constants for the land surface area. Although they are dependent on soil type, we believe that other land surface parameters will have a more pronounced effect on the hydrologic fluxes.
Continuous updating of g after each time step yields an updated spatial distribution of surface soil moisture. Note also that Equation (8b) can be used to determine the location of saturated areas. Saturated areas are those for which si ~< 0, or The growth and decay of saturated areas along stream channels can then be predicted by continuous soil moisture accounting, i.e. continuous updating of g. Continuous accounting of g is analogous to the accounting done in current GCM land surface parameterizations of grid average soil moisture.

a. THE SPATIAL DISTRIBUTION OF THE TOPOGRAPHY-SOILS INDEX AS A MODEL FRAMEWORK
The spatial distribution of the topography-soils index forms the foundation of our parameterization, through its control on soil moisture storage, water table depths, and its influence on near surface soil moisture. Infiltration and evapotranspiration capacities will be shown to depend on surface soil moisture, and thus the topography-soils index, in a later section. In this section, we discuss methods for quantifying the spatial variation in the topography-soils index. Once the variation in surface soil moisture has been quantified, the spatial variability in the water balance fluxes can be determined, and areal averaged fluxes can be computed. Dealing with the large quantities of topographic and soils data associated with large land areas can quickly become an overwhelming task, particularly when the global scale is considered. When the data is available, it is often in the form of a digital elevation model or a geographic information system (GIS). A common approach to utilize this data has been through the use of distributed models, where the land surface is treated as a collection of pixels whose size depends on the resolution of the data. The DEM represents the catchment digitally, and GIS data such as soil type and land use are layered over the DEM to represent land surface properties. The topography-soils index can then be calculated for each pixel. This approach works well at the catchment scale, but as larger scales are approached, computational and storage burdens can become excessive.
A simpler way to characterize the variability in the topography-soils index is by the use of probability density functions. In this paper we represent the large land area by its statistical distribution of the topography-soils index. The implicit assumption in this procedure is that knowledge of the exact pattern of the topography-soils index is not necessary to compute the land surface water balance fluxes. Research into a threshold scale, above which such exact patterns can be represented statistically, is in progress. This scale, called the representative elementary area (REA) by Wood et al. (1988), has been identified in the context of storm response modeling. Its existence in other branches of the water balance is currently being explored.
To perform water balance computations, the distribution is divided into a number of intervals. Each interval of the distribution represents the fraction of land surface area having that particular In(aTe/To tan t) value. A representative value of In(aTe/To tan/~) for the interval provides the necessary information for land surface hydrology calculations: for each interval value of In(aTe/To tan fl), the water table depth, soil moisture, and storage in the unsaturated zone are known. The corresponding land surface fluxes for each interval can be computed and weighted by the probability of occurrence of the interval. The areal averaged land surface hydrologic fluxes result. This averaging procedure is presented analytically in a later section as our land surface hydrology parameterization.
To draw an analogy to current land surface hydrologic 'bucket' models, this model will treat the grid square land surface as a distribution of buckets, whose capacity and surface soil moisture vary with topography, soil properties, and g. The resulting large scale hydrologic response is simply a weighted average of the responses of the individual buckets. Sivapalan et al. (1987) have fit a three-parameter gamma distribution to the In(a/tan t) values of two catchments in North Carolina. (Wolock et aL (1989) also fit gamma distributions to 145 catchments in the northeastern United States, and Wolock (personal communication), in other studies, has found that even the largest catchments fit a gamma distribution.) Assuming independence between transmissiv-ity and In(a/tan fl) values, Sivapalan et al. (1987) combine the gamma ~, ~b, X) for In(a/tan fl) with a normal (0, a 2) for ln(Te/To) and arrive at a gamma (/~*, 4~*, ;~*) for ln(aTe/Totanfl). They show that the parameters /~*, q~*, and X*, can be obtained from the parameters/~, q~, and ;( and an estimate of the variance of the saturated hydraulic conductivity. Therefore, only a topographic map or a DEM (to determine /~, ~b, and ;0, and an estimate of the variance of the conductivity are required to characterize the distribution of ln(aTe/T o tan p).
In the following sections we describe the appropriate equations for the water balance processes to be applied to each interval of the In(aTe/To tan fl) distribution. In a later section we present the analytical formulations of the averaging procedure which constitute the mathematical statement of our large scale parameterizations for grid square averaged land surface hydrologic fluxes.

a. INFILTRATION AT THE SOIL SURFACE
The governing equation for soil water flow in the unsaturated zone is given by Richards (1931) as where 0 is the moisture content, ~0 is the matric head, and K is the hydraulic conductivity. As is well known, the solution of (14) is not easy, due to the highly nonlinear nature of K(qJ), hysteresis, and the boundary conditions encountered in nature.  solved (14) with the simplified boundary conditions of an initially uniform moisture profile in the unsaturated zone, and a step change in soil moisture at the soil surface: His simplified infiltration equation is given by where f?~ is the infiltration capacity, S is the sorptivity (given by Sivapalan et al., 1987), Ks is the saturated hydraulic conductivity, and cKs includes the effect of graviy.
In this study, the spatial variation in infiltration rates arises because the sorptivity expression varies with ln(aTe/Totanfl) through its dependence on surface soil moisture. Consequently, ff will vary with each interval of the distribution. A mean areal value for Ks is used in calculating f~ for two reasons. First, the effect of spatially variable conductivity is already incorporated by the control of In(aTe/ To tan fl) on the local water balance. Second, we have assumed independence between hydraulic conductivity and the value of In(a/tan fl). For a given interval of the ln(aTe/T o tan fl) distribution, specifying a local value of Ks would imply knowledge of the covariation of hydraulic conductivity with In(a/tan fl). However, we have assumed that this covariance is zero. Figure la displays the variation of the depth of infiltrated water at capacity (Equation (15) integrated with At = 15 min) with soil moisture. Soil moisture is expressed as percent by volume, and the range of values shown extends roughly from the driest to the wettest conditions likely to be encountered on the land surface.
To couple the infiltration flux rate at the land surface to the atmosphere one must consider the rate at which rainfall is supplied to the soil. If the rainfall rate is very high, the infiltration capacity will be exceeded, and the actual infiltration rate will be equal to the capacity. However, if the precipitation rate is low, then all the rainfall can infiltrate. Therefore, the actual infiltration rate is given by where f,. is the actual infiltration rate, and P is the precipitation rate.
b. EVAPORATION FROM SOIL Two stages have been recognized in the unsteady drying of a soil profile (see Brutsaert, 1982, andHillel, 1980). In the first stage, the moist soil profile has no problem supplying all the water that the atmosphere demands. Thus, this stage is known as the atmosphere controlled stage. Evaporation proceeds at the potential rate, which is dictated by external climatic conditions. The duration of this stage depends on the rate of the atmospheric demand and the ability of the soil to supply moisture at this rate. Hillel (1980) notes that this stage is frequently brief, and usually ceases within a few days. As the soil near the surface dries out, moisture can no longer be delivered at the rate demanded by the atmosphere. Instead, the moisture delivery rate is limited by the properties of the soil profile. Thus, this stage of soil drying is known as the soil controlled, or falling rate stage. Brutsaert (1982) notes that at any one point, the transition from soil to atmosphere control is rapid, but over the entire catchment, the changeover will be gradual.
The governing equation for the soil controlled stage of evaporation is obtained by combining soil water continuity with Darcy's law, which yields A simplified formulation considers the soil controlled stage as a desorption problem only. Neglecting gravity, (17) For the simplified boundary conditions where 0 d is the moisture at the dry soil surface, (18) can be solved using the Boltzman transform (~b = zt 1/2). Equation (18) reduces to an ordinary differential equation, and the evaporation capacity, j'~, is given by where D is the desorptivity (Sivapalan, unpublished), which is dependent on soil type, Od, and the surface soil moisture, 0;. In Equation (19), the desorptivity varies with soil moisture content and thus In(aTe/To tan fl). As in the infiltration case, will vary with each interval of the In(aTe/To tan fl) distribution. Figure lb displays the variation of the depth of evaporated water at capacity (Equation (19) integrated with At = 15 min) with soil moisture. As in Figure la, the range of soil moisture values shown extends roughly from the driest to the wettest conditions likely to be encountered on the land surface. Coupling this evaporation expression to the atmosphere requires a knowledge of the atmospheric demand for water vapor. When the atmospheric demand, or potential evaporation, can be met by the local land surface area, then the actual evaporation rate is equal to the potential evaporation, and that fraction of the land surface area is subject to atmosphere controlled evaporation. When the local land surface area can no longer meet the atmospheric demand for water vapor, the actual evaporation rate is equal to the evaporation capacity, and that fraction of land surface experiences soil controlled evaporation. At any time, the actual evaporation rate can be expressed locally as fe = min [f~e, epe ] (20) where fe is the actual rate of evaporation, and epe is the potential rate. The potential evaporation is assumed known from atmospheric data or GCM variables.

C. TRANSPIRATION BY VEGETATION
The importance of vegetation in the hydrologic cycle cannot be underestimated. Molz (1981) points out that well over half of the water returned to the atmosphere as evapotranspiration flows through the soil-plant system. It has also been recognized that plants are not simply passive wicks conducting moisture from soil water reservoirs to the atmosphere. Rather vegetation is a dynamic component of the soil-plant-atmosphere system, which actively regulates its internal mechanisms (e.g. stomatal closure) in response to changing soil and atmosphere conditions. The current state of knowledge of transpiration processes should be reflected in any model of the vegetative component of the hydrologic system. As with the storm  response model, our approach is to keep the transpiration component simple, while maintaining the appropriate level of the process physics, which are briefly described below. Denmead and Shaw (1962) are credited with first confirming the effects of dynamic soil and atmosphere conditions on transpiration rates. They showed that the ratio of actual to potential transpiration depended on both the potential rate of transpiration and soil moisture content. In dry soils, vegetation can transpire at the potential rate if that potential is low; in wet soils, actual transpiration can fall below the potential rate if that demand is too high.
We have chosen to model the uptake of soil moisture by roots using the extraction function of Feddes et al. (1976). A macroscopic approach is assumed, where the entire root system is viewed as a diffuse moisture sink, rather than considering flow to individual roots in a geometrically complex root system (microscopic approach). This extraction function is shown below. S is the actual rate of transpiration, and Smax is the maximum rate at which the vegetation can deliver moisture to the atmosphere. Smax is taken as the potential rate of transpiration, ept, and is assumed known from atmospheric data or GCM variables. The value of surface soil moisture at which transpiration can no longer be sustained at the potential rate is expressed as 01, and 02 is the wilting point for vegetation. When 0i ~>01, transpiration occurs at the potential rate. When 02 ~< 0i < 01, vegetation transpires at a maximum sustainable rate, or capacity, which depends on the availability of soil moisture. When 0i remains below 02 for extended periods, wilting follows. The supply and demand dynamics described by Denmead and Shaw are incorporated by allowing the value of 01 to vary with atmospheric demand and surface soil moisture, so that at low levels of potential transpiration, even dry areas can transpire at the potential rate, etc. For simplicity, this value of 01 is taken to be the same value at which the change from atmosphere controlled to soil controlled evaporation occurs. The wilting value of soil moisture, 02, is however taken as a constant.
For each interval of the In(aTe/To tan fl) distribution, 01, the surface soil moisture, will vary. Consequently, during any time step, a proportion of the vegetated land surface will transpire at the potential rate, and a proportion will transpire at the transpiration capacity, i.e., a fraction of the vegetated land surface will be subject to atmosphere controlled transpiration, and another fraction will be under 'vegetation controlled' transpiration. In areas where soil moisture conditions are exceedingly dry, an additional fraction of vegetated land surface will experience wilting conditions. Figure

Saturation Excess Runoff
Given a value of ln(aTe/T o tan fl) for a particular interval of the distribution, the available soil moisture storage can be calculated using (8). Rain falling on saturated areas (si ~< 0) is transformed immediately into runoff. Additionally, over parts of the land surface where the available soil moisture storage is small, for example, close to streams where the water table is high, the local storage deficit may be satisfied during a storm. To determine when saturation excess runoff occurs on these areas, one need only keep track of the volume of infiltration and compare it to s;. When the infiltrated volume of moisture exceeds available soil moisture storage, saturation excess runoff occurs on that interval of the distribution.

Infiltration Excess Runoff
Infiltration excess runoff occurs when the local infiltration capacity is exceeded by the precipitation intensity, or P >3~. That rainfall that is not infiltrated runs off downslope to the stream channel.

Subsurface Flow
The contribution of the land surface to runoff by subsurface flow, Qs is given by integrating the downslope saturated flow, q; along both sides of the stream network: Os = fq; dL (22) JL or where Qo = ATe exp( -2) (24) (Sivapalan et al., 1987) L is the total length of the channel network (both sides), and A is the area of the land surface. Again, this is the subsurface flow exiting the hillslope at the stream channel. Ncte that the quasi-steady state approach implicitly incorporates the dynamics of subsurface flow on hillslopes during storms; as g is updated, the water table profile shifts in response. The parameters m and Qo are physically based, and can be obtained from field and map information (Beven and Kirkby, 1979;Beven et al., 1984) Alternatively, they can be determined by calibration to a number of recession curves (Beven and Wood, 1983). The average storage deficit at the start of a simuiation can then be obtained by inverting (17).
For large land areas and particularly for GCM grid squares, appropriate streamflow data is most likely not available. Research is in progress in which we attempt to extract these parameters from remotely sensed soil moisture.

The Spatial Distribution of Precipitation Intensity
Given the large areas under consideration, subgrid scale rainfall intensities will vary considerably. When using spatially distributed models, this poses no problem: rainfall inputs can vary with each node or pixel in the model. Space-time varying inputs can be determined from radar estimates of rainfall, stochastic rainfall models, or from weighted point measurements from raingauges. However, in a lumped modeling approach such as this, only one value of rainfall can be used as input. Hydrologists realize that such an average estimate of rainfall cannot possibly capture the natural spatial variability of the rainfall process. The experience of climatologists suggests that the one value of precipitation computed for a GCM grid square is unrealistically low.
One way to reconcile this disparity is to assume that the rainfall intensity predicted by the GCM, or averaged from ground data, is the mean of a distribution of possible intensities falling on any one location. We assume that the point rainfall intensity, P, is exponentially distributed (Eagleson et al., 1987), with mean E[P], as In this manner, each interval of the ln(aT~/To tan r) distribution, i.e. each fraction of the land surface, can be subjected to a range of possible rainfall intensities, each with an associated probability of occurrence. For each interval of the In(aTe/ To tan r) distribution, an expected response to the distribution of rainfall intensities will be calculated. These interval responses will themselves be aggregated to yield the expected or areal averaged storm response for the land surface.

Macroscale Parameterizations for Land Surface Water Balance Fluxes
In this section we domonstrate how to aggregate the water balance fluxes of the individual intervals of the In(aTe~To tan r) distribution. This procedure is equivalent to calculating the areal averaged fluxes for the land surface. Computationally, fluxes for each interval are weighted by the corresponding probability of the interval, and large scale average fluxes result. The analytical formulations are presented here. These constitute our land surface parameterizations for storm runoff and evaportranspiration processes. The equations are proposed for use in GCMs or for use in a large scale water balance model. Calculations require minimal computer time, and only a knowledge of the topography and some soils and vegetation properties is needed in addition to routine land surface information. In the equations presented, the distributions of precipitation and the topography soilsindex are assumed to be independent. Also, flux rates for precipitation, infiltration, evapotranspiration, and potential evapotranspiration, are assumed to now represent a depth over some specific time interval.
Equation (13) gives the relationship between the topography-soils index and saturated areas: those points with a topography-soils index greater than or equal to the right side of (13) are saturated. For convenience, let us define a 'saturation' value of the topography-soils index, ln*(aTe/T o tan/~). This is the value at which the storage deficit, si, is just equal to zero. Topography-soils indices greater than or

I. The Expected Value of Infiltration Excess Runoff
The expected value of the depth of infiltration excess runoff for the large land area, where x = ln(aTe/T o tan fl). Equation (27) states that on unsaturated areas of the land surface (0 ~< x < ln*(A Te/To tan fl)), infiltration excess runoff occurs when the local infiltration capacity has been exceeded. As previously mentioned, infiltration capacity is parameterized here as a function of surface soil moisture which depends only on the local value of the topography-soils index.

The Expected Value of Saturation Excess Runoff
All rain falling on saturated areas of the land surface is transformed immediately into saturation excess runoff. Additionally, over some fractions of the land surface, storage deficits will be satisfied during a time step. Rain falling on these newly saturated areas will also be transformed into saturation excess runoff. j ~X = ln*(aTE/T 0 tan fl) /~ = oo + Je (P -Si)fp(P)fx(x) dp dx. X=0

=S i
The first term on the right hand side represents runoff generated on those areas that are saturated at the start of a time step, while the second term represents the runoff generated on those areas that become saturated during a time step.

The Expected Value of Total Storm Runoff
The expected value of the depth of total storm runoff, E[Qt] is given by summing (27) (29) for any time step.

b. EVAPOTRANSPIRATION
The long term distribution of vegetation type and density depends on the interaction of climate, topography, and soils. We can observe, for example, that vegetation density increases along stream networks where the water table is high (high values of In(aTe/To tan//)). We can hypothesize that vegetation type, and thus wilting point, 02, will vary locally with ln(aT~/To tan fl), and regionally with the large scale average soil moisture deficit, ~. At present, the spatial variation in these parameters can be quantified somewhat by remote sensing. Nevertheless, their relationship to the distribution of the topography-soils index has received little attention. Extracting this information from remotely sensed data and investigating these relationships is an area worthy of future research: however, for the purposes of this paper, some simplifying assumptions will be made. First, we assume that each interval of the In(aTe/To tan//) distribution is covered by a constant vegetated fraction, F,. Second, as previously mentioned, we assume that one value of 02 applies to the entire vegetated canopy. Now we can proceed with our land surface evapotranspiration parameterizations.

I. Areal Averaged Bare Soil Evaporation
Saturated areas will evaporate at the potential rate. The areal averaged depth of bare soil evaporation by saturated areas is given by where the subscript s implies saturation and the subscript ac implies atmosphere control.
In addition to saturated areas, at each time step there will be areas on the land surface that are unsaturated, but can still supply moisture at the rate demanded by the atmosphere. The contribution of these areas is given by where the subscript u implies unsaturated areas, and ln'~(aT~/To tan//) is obtained by setting epe=~ solving for Oi, the surface soil moisture (which is a parameter of D) and inverting (12), (5), and (8). Because the parameterization of D employed is unpublished, we cannot present an expression for lntS(aTe/To tan fl). Still, we point out that this value of the topography-soils index is effectively a threshold between atmosphere and soil controlled areas of bare soil evaporation (the superscript ts indicates a threshold for soil control.). The threshold value depends on both soil moisture conditions and the atmospheric demand, so that for any time step, the total average depth of evaporation contributed by atmosphere controlled areas is given by = lntS(aTe/T 0 tan B) The areal average depth of evaporation contributed by soil controlled areas is given by

Areal Averaged Transpiration by Vegetation
All saturated areas will transpire at the rate demanded by the atmosphere. The areal average depth of transpiration by these areas is given by where again the subscript s implies saturation and the subscript ac implies atmosphere control. Analogous to the bare soil case, there will be land surface areas which are not saturated, but can still maintain transpiration at the potential rate. The threshold topography-soils index distinguishing 'vegetation' controlled transpiration from atmosphere controlled transpiration lntV(aTe/To tan fl) is for simplicity assumed equal to the threshold index for the bare soil case, lnt~(aTe/To tan fl). Those areas that are unsaturated and transpiring at the potential rate yield an average depth of T.. = at= ln'~'(.r~/ro t~. ~) ePOc~':(x) dx (36) where again the subscript u implies unsaturated areas.
Unsaturated areas that are unable to meet the atmospheric transpirational demands, i.e. those land surface areas under plant controlled transpiration, produce an average depth given by ~x = lntV(aT e/T O tan ,g) T,~c = Jx=l nW( aTe/ Tot an~) | Sfx(X) dx (37) where the subscript vc implies vegetal control, and lnW(aTe/T O tan fl) is the value of the topography-soils index where 0j is equal to the wilting value, 02, and is obtained from Equations (12), (5) and (8).
Total transpiration for the vegetated land surface area, T,, is given by This expression is obtained by summing (35), (36), and (37).

Areal Averaged Evapotranspiration
Finally, areal average evapotranspiration, ET is given by

Model Operation
The first step in model operation is to estimate the areal averaged initial soil moisture deficit, g. This can be accomplished by inverting (23), or with the help of remotely sensed soil moisture, or from measurements (since g is a function of average water table depth by (5), regional estimates of average water table depth would suffice). Since g implies a distribution of surface soil moisture through Equations (8), (5), and (12), then the areal averaged fluxes of runoff and evapotranspiration can be computed, given the atmospheric forcing and Equations (29) and (39). These fluxes can then be used to update L and the procedure is repeated for the next time step.

r RUNOFF RATIO AND EVAPORATION AND TRANSPIRATION EFFICIENCY
Manipulations of the equations presented above yield the dimensionless runoff ratio and evaporation and transpiration efficiencies. The ratio of surface runoff to precipitation, or the runoff ratio, R = E[Q,]/E[P] is given as

R
where As/A represents the fraction of saturated land surface. In this expression, we ignore the contribution of subsurface flow to total runoff, and consider only the ratio of surface runoff to precipitation.
The bare soil evaporation efficiency, fls = Eb~/epe, given by dividing Equation where A t ..... /A represents that fraction of land surface not saturated but still able to evaporate at the potential rate. The subscript 'trans' implies a 'transitional' region where soil moisture decreases from saturation values to a threshold value below which evaporation proceeds at the soil controlled rate. The subscript s represents bare soil.
The transpiration efficiency, fly = Tt/ep, is obtained by dividing Equation (38)

ep,
Analogous to the bare soil case, A t .... ,,/A represents that fraction of land surface area not saturated but still able to transpire at the potential rate, the subscript 'trans' implies 'transitional' between atmosphere controlled and vegetation controlled transpiration, and the subscript v represents vegetation.

Discussion
The runoff ratio and evaporation and transpiration efficiencies effectively characterize the large scale hydrologic response. The various modes of land surface behavior in response to atmospheric forcing are represented by the individual terms in each of the ratios. The evaporation efficiency consists of a term representing the fraction of saturated land surface, where fls is 1, a term representing a transitional region where the soil is unsaturated, but the evaporation efficiency is still 1, and a soil controlled area term. The transpiration efficiency is composed of analogous terms to fls : a saturated area term where transpiration is at the potential rate, a transitional area term where transpiration is still at the potential rate but the soil moisture is decreasing, and a vegetation controlled term. The runoff ratio consists of a saturated area term, where runoff is a maximum, and although a transitional area term does not fall out neatly from the parameterization, it is in fact contained within the infiltration excess runoff term in (40). It is easy to conceptualize that as surface soil moisture decreases from saturation values to drier values, a threshold value of soil moisture will be encountered, beyond which infiltration capacities are too high to produce runoff. Over the range of soil moisture between saturation and this threshold value, runoff will decrease from a maximum rate to zero. The fraction of land surface with surface soil moisture lower than this threshold value will produce no runoff. Because E [Qinf] in (28) is an integral over all unsaturated areas, both the transitional region term and the drier region term are contained in the infiltration excess runoff term in (40).
In Figure 2a-c the runoff ratio, R, is plotted versus g, the average storage deficit, for increasing levels of rainfall intensity. The range in magnitude of g can be considered due to seasonal variation, with higher numbers representing drier conditions. The solid line in each plot represents the runoff ratio. For one value of the average storage deficit, R increases as the rainfall intensity increases from Figure  2a to 2c. There are two explanations for this. First, infiltration excess runoff is occurring over a large fraction of the land surface as the rainfall intensity increases. This is shown by the difference between the solid line (R) and the dashed line (A s/A)  explained by the fact that less of the land surface can supply moisture to the atmosphere at the potential rate as that potential increases. This fact is demonstrated by studying the middle line in Figures 3a through 3c for one value of g. This line is the sum of the first two terms on the right hand side of (41), i.e. this line represents the total fraction of land surface under atmosphere controlled evaporation. As the atmospheric demand increases, less of the land surface can supply moisture at the potential rate, and at the highest levels of demand, only saturated areas and another very small fraction of land surface area can contribute moisture at the rate demanded by the atmosphere (Figure 3c). Similar results are shown in Figures 4a-c for the transpiration efficiency, although for the same level of atmospheric forcing fly is greater than fls. In any of Figures 3a-c, as g increases, fls generally decreases. As soil moisture conditions dry out, less of the land surface is able to supply moisture to the atmosphere at the potential rate. Again this is shown by the middle line in any of Figures 3a-c. As g increases, the proportion of land surface area under atmosphere control decreases. The same trend is evident in Figures 4a-c for the transpiration efficiency.
The parameters used to generate Figures 2-4 are given in Table I.

Summary
A macroscale model for the land surface hydrologic fluxes is presented. The model is proposed for use in atmospheric GCMs to improve upon the current simplified land surface hydrology parameterizations. Additionally, the model functions independently of a GCM as a large scale water balance model. The model incorporates subgrid scale spatial variability in topography, soils, vegetation, and rainfall to predict the space-time distribution of surface soil moisture over a large catchment or GCM grid square. Given the distribution of surface soil moisture and the associated probabilities, the distribution of surface runoff, evaporation, and transpiration rates can be computed, and the areal averaged hydrologic fluxes can be determined by integration over the distribution of soil moisture.
The parameterization presented here suggests that the interactions between the land and the atmosphere that produce the macroscale water balance fluxes are dominated by three broad divisions of the surface soil moisture distribution: saturated, transitional, and relatively dry. The relative magnitude of the saturated area varies seasonally, and on the shorter time scale of storms. The magnitudes of the transitional and relatively dry areas is a function of both climatic forcing and the state of the soil surface, and the threshold values of soil moisture between regions will vary with flux type. Saturated areas contribute runoff and evapotranspiration at the maximum rate. Transitional areas contribute runoff at a decreasing rate (spatially) as the soil moisture decreases (and infiltration capacities increase) within the boundaries of this region. Transitional areas contribute evapotranspiration at the maximal rate. Relatively dry areas contribute no runoff and contribute evapotranspiration at a decreasing rate (spatially) as the soil moisture decreases within these regions.
Embedded in this parameterization is the subgrid scale variability in hydrologic processes and land surface properties which we believe are crucial to the dynamics of land surface/atmosphere interactions. In response to these water balance dynamics, the saturated and threshold areas will expand or contract diurnally and seasonally. It is precisely this subgrid scale heterogeneity in surface properties and dynamics that we believe are important in determining the large scale averaged response. The temporal variation in the terms of the above ratios, their sensitivity to various climates, vegetation, and topography-soils characteristics, and validation of the model on large watersheds, are the subjects of future research.