Sierra Nevada Ecosystem Project: Final report to Congress, vol. III. Assessments, Commissioned Reports, and Background Information. Davis: University of California, Centers for Water and Wildland Resources, 1996.  Author for Correspondence:   Susan L. Ustin   Department of Land, Air, and Water Resources   University of California   Davis, CA 95616   Phone:  (530) 752-0621   FAX:  (530) 752-5262   email:  slustin@ucdavis.edu

INTRODUCTION

Hydrological processes are influenced by factors like climate, topography, soil type and structure, and vegetation. Because these factors vary continuously in the landscape and interact together, it is difficult to predict the impacts of changes in any of their properties within the watershed without a simulating tool. The Surface-Subsurface Model is coupled with a Geographic Information System (GIS) that provides an easy way to manipulate large arrays of distributed data. Thus by changing the inputs and running the program, we can visualize and analyze their effects on the outputs, spatially and temporally. Finally, this model can be considered a management tool to predict the impact of landscape transformations like forest fires, logging, grazing or other activities on the hydrologic budget. The model can simulate the effects of complex interactive phenomena. It may improve our understanding of the hydrologic processes in a mountainous area.

What effect will changes in the water regime have on terrestrial ecosystems? This is a central question affecting predictions of ecosystem response to global climate change. Understanding what components of the hydrologic cycle and related ecosystem processes are sensitive to land use and climate change is critically important for the entire Sierra Nevada ecosystem. A far better understanding of regional hydrologic processes is needed if we are to understand how land management interacts with climate so that realistic strategies can be developed for balancing the demand for urban and agricultural water use. The interaction between climate and hydrologic systems is through plant water use. The biology of terrestrial and aquatic ecosystems depends on an integration of the quantity, quality, and temporal delivery of water resources, the complexity of which requires an integrated analysis of surface and subsurface processes and conditions. The distribution of terrestrial ecosystems is primarily determined by the interaction of climate and topography and the type and frequency of the disturbance regime on the landscape. The extreme Mediterranean climate in California and the frequency of fire create a highly dynamic landscape, where relatively small climate shifts can have significant impact on the local distributions of biome types. Water uptake and transpiration rates determine the amount of CO₂ available for photosynthesis, the mineral nutrition status of the plant, and limits the potential growth rate. Soil moisture, nutritional status, and disturbance regime (type and frequency) largely determine the impact of environmental "stresses." Population dynamics interfaces with the response to environmental stresses and determines whether directional succession will lead to changes in the distribution and abundance of plant communities on the watershed.

To a first approximation, canopy water use is driven by the surface energy budget, available soil moisture and the foliar biomass or surface area. The seasonality and distribution of soil moisture determine both the type of vegetation and the biomass accumulation. Several years of below normal precipitation can cause increased mortality to the forest directly, or indirectly, by increasing susceptibility to insects, pathogens, and wildfire as occurred during the drought of 1987-1993. If demand for latent heat exchange changes (due to a change in net radiation, either total or seasonal), or the supply of water (soil moisture either total or seasonal), then species distributions and/or densities and biomass will change. It is possible to predict the direction of changes in vegetation type or density under different climate change scenarios provided sufficient knowledge of seasonal soil moisture demand and energy budgets by different community types are understood. The ecosystem implications of climate change can be predicted through changes in evaporative demand, phenology, life history, and ultimately shifts in species distributions and biome boundaries. Such changes could be predicted by the ecohydrologic model.

CAMP CREEK CASE STUDY

A smaller sub-watershed of a tributary into Camp Creek was used for hydrologic simulations and scenarios. The Snow Creek and Pebble Creek watersheds were selected by the SNEP team for several case study scenarios. Snow Creek is centered at 38° 42' latitude and 120° 25' longitude and is 3.55 km². Pebble Creek is centered at 38° 42' latitude and 120° 26' longitude and is 2.81 km².

Another advantage in studying Camp Creek is the availability of satellite data to drive the hydrology model. Satellite data provides a basis for enriching the spatial information in the model by making it possible to develop a truly spatially distributed hydrology model. The hydrology model is now running at the 30 meters spatial scale of TM images on the Camp Creek drainage. There is an excellent GIS database for this watershed that includes Landsat Thematic Mapper satellite images, DEMS, digital soil and geology maps, roads, forest vegetation and inventory maps, county build-out maps, and many other data layers in the SNEP database. Data from U.S. NOAA Climatological Weather Station data, California Climate Data Center (CDEC), California Irrigation Management and Information System (CIMIS) data, and data are from the forest service are also available. Quarter hour water level recording data was available for intermittent periods between January 1987 and June 1990 at three stations (Darlington Canyon, Pebble Canyon, and Little Light Canyon) within the watershed. Because there is no recorded stream cross section the data can't be converted to a hydrograph. We searched for existing hydrologic and weather data which could be used for model validation. Well data from four stations in Darlington Creek (Section 26, TIOR14) were available from the Eldorado National Forest which were installed by the forest hydrologist and geologist to monitor herbicide transport with sub-surface flow. The wells are 38, 22, 65 and 16 ft deep, respectively. Well logs were available with water depth recorded manually beginning in Sept., 1991 and continued weekly through Sept. 29, 1992. Continuous recording rain and stream gauges were installed to complete the hardware for monitoring the hydrologic cycle and battery operated data loggers were installed. However, because of intermittent power failures, only incomplete data were provided. The forest also recorded tipping bucket rain gauge data at about 2780 in elevation but which was unreliable because of uncertainties during snow fall periods. Due to the limited nature of the hydrologic data for Darlington Canyon we excluded it from use in the current validation of the Ecohydrologic model. However, continued investment in the site would be useful for future validation work and the data stations should be maintained. Using recording devices designed for a wider range of weather contions may be an important upgrade for future studies.

The Eldorado National Forest also has adequate field-measured ecologic unit site data within and around the watershed to evaluate data quality for soils, geology, potential natural vegetation, and special features. The Forest Inventory Analysis (FIA) for Eldorado National Forest was used to evaluate vegetation characteristics and was a database on 26,000 trees from 241 plots.

The general structure of the model is shown in Figure 2. This schematic figure of a watershed illustrates the spatial relationships among the cells in the model and shows the connection between the inputs and outputs for a cell. The model was run with cells having a size of 30m by 30m for the Camp Creek watershed. Another way to visualize these relationships is shown in Figure 3 which illustrates the inputs and outputs for a cell in the watershed. This figure describes some scenarios that might be considered using the model, such as climate processes, wildfires, and forest management and illustrates some environmental effects that might be evaluated.

OBJECTIVES

TERRESTRIAL ECOSYSTEMS AND HYDROLOGIC MODELS

It is necessary to predict spatially distributed patterns to predict the impact of possible climate variation, particularly if the impact associated with global warming is to be understood. It is difficult to predict climate changes such as those that might result from atmospheric trace gas forcing. Long-term GCM estimates suggest that one effect of doubling the present atmospheric C0₂ concentration could be wan-ner winter minimum air temperatures, which would lift the snow line several hundred meters. Predicting the possible distribution of snowpack under these scenarios is difficult but important if hydrologic response is to be evaluated. These warmer winter storms would lead to more precipitation falling as rain and less as snow even if the total amount of precipitation did not change. Changing the distribution of precipitation would increase the probability of winter floods and decrease montane snow storage, the hydrologic component that is critical for the late summer water. supply. Either climate outcome would bring adverse effects to terrestrial and aquatic systems in the Sierra Nevada. Understanding hydrologic process and the hydrologic response to land use and climate change requires an integration of all hydrologic components. It is not enough to develop the most inclusive and theoretically complex hydrologic model, because the input parameters need to be easy to derive from existing data sources and realistic. So some balance is needed between rigorous physically driven processes, scale of the processes, and data availability.

Most "physically based" hydrology models fall short in combining surface and subsurface flow much less including vegetation effects on the water budget. The technical meaning of physically based models have various interpretations. Generally what is meant is a "quasi-physical" model that uses simple physical equations for hydrologic phenomena. It must be kept in mind that the physical laws behind such models have been validated under specific and narrowly defined criteria and the linearity of the processes and appropriateness of multiple time/space scales may not be valid when applied over entire watersheds. Nonetheless, these models have largely been focused on regional scales while resources managers typically plan at a watershed scale. Furthermore, models do not fully exploit remote sensing data--the only near-contemporaneous spatially detailed data that exists.

Surface water simulators such as US Army Corps of Engineers TABS-MD calculate surface flow to streams but ignore the groundwater system. Similarly, groundwater analysis models (e.g. the USGS MODFLOW) assume given or static inputs from the surface as boundary conditions which disconnected them from the vadose zone and vegetation effects. Jordan (I 992) applied two of the best developed models, the SHE (Systeme Hydrologique Europeen; Abbott et al., 1986) and TOPMODEL (Topographic Models; Beven and Kirby, 1979), but could not model fast subsurface flows despite including all physical processes in the model. The SHE model used cells of 250m by 500m (1% of watershed). Gao et al. (1993), developed a fine scale (15 m by 15 m) spatially explicit physically-based surface-subsurface model to overcome the limitations of surface only or subsurface only approaches. They used a 1D Richard's equation to simulate unsaturated flow, a 2D groundwater flow equation assuming horizontal flow, a 2D kinematics wave equation to simulate overland flow, and stream flow was routed using the St. Venant equation. It was possible to model the basin solving the partial differential equations (PDES) implicitly because (1) the watershed was only 4.4 ha. and well-instrumented, (2) the groundwater flow system was deep and isolated from the surface and as such was ignored, and (3) stream flow was not modeled. Even with these simplifying assumptions, the complexity of their approach is impractical at the scale of typical watersheds nor is the essential data to derive the model available. The model is computationally intensive due to the small time steps and very large matrix inversions that are necessary to solve implicit solutions of PDEs over large areas using grid cells of realistic size. Furthermore, the physically based models (described above and others) input parameters are largely depend on field observational data sources and existing maps. Finally, the interactions between hydrologic properties and ecological conditions (e.g., evapotranspiration rates and patterns) are not examined.

Even fewer models consider interactions of the physical hydrologic system with ecological conditions, despite obvious connections between plants and the interception of precipitation, soil water retention, and evapotranspiration. Hydrological models such as Regional Hydroecological Simulation system (RHESSys) are somewhat more biologically based but lack physical rigor. In recent years several models that address spatial and temporal information have been developed, such as the SiB2 model (Sellers et al., 1986), MTCLIM model (Running et al., 1987; Hungerford et al., 1989), and Forest-BGC model (Running and Coughlan, 1988), that can simulate climate and some hydrologic components for large scale studies. The minimum spatial scale in these models is on the order of kilometers and the temporal scale on the order of days. A distributed hydrology-vegetation model was developed for complex terrain (Wigmosta et al., 1994) that used 180 m spatial and one-day temporal scales to study surface and subsurface water movement. However, these models were developed for ecological studies so the main focus has been to predict evaporation and plant transpiration while ignoring both infiltration and surface flow. A simplifying assumption has been that the study area had a homogeneous soil and the surface cover parameterization used a preexisting vegetation map that may not correctly represent current vegetation distributions. The objectives, input parameters, and spatial and temporal resolution of these models have limited their application when applied to watershed scale hydrologic studies.

ECOHYDROLOGIC MODEL

MODEL DESCRIPTION

Within each cell, all surface and subsurface flows are linked and water flows laterally between cells. Surface flow is stimulated as steady during a time step using Manning's equation. The equation is applied in the x and y directions of the horizontal plane and is a simplification of the two dimensional solution of the St. Venant equations. Forms of this equation have been widely applied in surface flow models (Abbott et al., 1986; Jordan, 1992; Gao et al., 1993). Surface roughness is derived from remote sensing data and available vegetation maps. Surface slope is derived from DEM data, and the water ponding depth is estimated from continuity of water movement. Infiltration is calculated using the method developed by Green and Ampt (1911). The parameters are all derived from available soil survey databases (e.g., USDA SCS or USDA Forest Service) and derived relations using standard equations (Maidment, 1993). A one-layered soil profile (surface undecomposed organic matter, root zone/unsaturated zone, and saturated zone) is used for calculating soil moisture distribution, unsaturated flow, and saturated subsurface flow. The horizontal water movement in the subsurface and ground water recharge are calculated from Darcy's Law. In the calculation for subsurface flow and determination of water table, the Dupuit assumption is assumed to be satisfied.

The study area is gridded at the same scale in X and Y directions (Figure 2). We assume the land surface and subsurface physical properties are homogenous within each grid cell even though the properties are heterogeneous over the whole study area. Figures 2 and 3 illustrated how the grid and the model are abstracted from the natural landscape. A mass balance is computed for each cell in the grid. Figure 5 shows a three dimensional representation of the hydrologic processes in the model and an X- Y plane view of the surface and subsurface flow paths. This figure describes the inputs and outputs in terms of the hydrologic fluxes.

The spatial scale of the model can be varied from meters to kilometers (depending on input data resolution) and temporal scales from seconds for storm studies to daily for ET and subsurface flow. This model is suitable for application to any study area because the model is not trained on empirical data. Because this model is physics-based, and uses remotely sensed data in a GIS database, it can be applied to any study area and produces surface and subsurface hydrographs as well as soil moisture history for all cells in the problem domain. All of the parameters and variables needed by this model are directly or indirectly found from satellite remote sensing data, USGS topographic data, available soil survey data, and meteorological data. Land use and construction layers (forest, roads, stream locations, buildings, and free surface water) are classified from Landsat TM data while others are from SNEP databases.
 

Climate Variables.

Precipitation

NOAA Climatological Data from sites identified in Table 1, were used to compute linear and third order monthly regressions for precipitation rates. We avoided stations with shorter data records and those located in sites where rainfall rates were possibly biased by their location. Since rainfall is strongly dependent on topographic elevation, we correlated rainfall and elevation by developing an empirical "environmental precipitation lapse rate" for the Camp Creek region. Because the weather stations are located approximately along the same latitudinal gradient, we assume that latitude is implicit in the function. In an initial examination of these data, we observed that rainfall generally increased with elevation up to about 2000 m and declined at higher elevations. Because of this nonlinearity, we computed two different monthly correlations, linear and third order (Table 2). The minimum monthly correlation coefficient for predicting rainfall as a function of elevation, for a third order regression was 0.86 (October). The best-fit 30 year mean monthly precipitation predictions for February (Figure 6a) and May (Figure 6b) are shown as a function of elevation.

The fit of the data set were tested using a separate data set, a recent eight-year record (1987-1994) with both hourly and daily data obtained from 16 NOAA Climatological Data stations and the NCAR climate database. These sites crossed the same elevation gradient (Table 3) and latitude as the long-term dataset used to create the relationships. These data also show the that peak rainfall occurs in the mid-elevation range around the 2000 m contour. This result is consistent with recent model predictions for orographic precipitation (1994). Thus the short-term record from 16 stations located near the watershed displayed monthly precipitation and elevation patterns that were similar to the long-term mean precipitation patterns. We did not use these data directly because the short-term record was more noisy than the 30 year pattern. Because the short-term record spanned the drought years they were anomalous in their more sporadic rainfall distributions. A second reason for the greater noise in the eight year data set is because data records from several stations were incomplete and the missing data could also have biased the results. Again, first and third order regressions for monthly means showed good correlations, with regression coefficients greater than 0.87.

The final data set used to evaluate weather generation is the California Irrigation Management Information System (CIMIS) weather station at Camino, CA and located within the Camp Creek study area. The regression between predicted and measured cumulative annual rainfall for the 1995 water year are shown in Figure 7. This figure provides an independent assessment of the accuracy of the weather predictions for a random location in the watershed and illustrates the reliability of the climate sub-model. The first order regression prediction closely approximated the measured precipitation at the Camino station. The third order regression prediction was slightly less good in this test of the prediction.
 

Net Radiation

Net shortwave radiation for each cell was simulated in the model (following the methodology in MT-CLIM, Running and Coughlan, 1988) and was used to derive daily radiation budgets. These measures were adjusted for cloudiness using the Campbell and Bristow method and extrapolated spatially based on the solar zenith track for each time of day and date and the topographic conditions in each cell. Daily net shortwave radiation was modified by calculating an east and west horizon mask.
 

Energy Budget Equation

Vapor Pressure Deficit

Vapor Pressure (kPa) under ambient conditions is
 

(2)

where

T_site =mean daylight site Temperature ° C

The saturated vapor pressure in the atmosphere is determined by
 

(3)

where

e_a is the vapor pressure in the atmosphere
e_s is the saturated vapor pressure in the atmosphere at that temperature
T_dew = dew point, ° C

The vapor pressure deficit VPD (kPa) is calculated as
 

VPD ea – es

(4)

 

Albedo

Snowfall

Snowmelt

We have incorporated an energy-balance method that estimates radiative transfer and heat advection to the snow, and assumes that the other fluxes are relatively small and can be neglected. Radiative transfer is the most important flux driving snowmelt in the Sierra Nevada (Kattelmann, 1990). Estimation of heat advected to the snowpack is important for prediction of rapid snowmelt due to rain-on-snow storm events.

Radiative transfer calculations include daily short- and long-wave radiation exchange with the snowpack. The snowpack is always assumed to be below the canopy. Net shortwave radiation at the top of the canopy is predicted as described above (Net Radiation). Incoming long-wave radiation flux is a function of atmospheric emissivity, air temperature, and estimated cloudiness.

Net shortwave radiation below the canopy (transmitted) is calculated by attenuating above-canopy shortwave radiation (incident) as a negative exponential function of leaf area index. If leaf area index is less than 1.0 then all incident radiation is assumed to be transmitted. Below canopy radiation is calculated as
 

It = I0 e-k*LAI

(8)

where I_t is transmitted shortwave radiation, I₀ is incident (above-canopy) shortwave radiation, LAI is one-sided leaf area index and k is an extinction coefficient (0.5).

Incoming longwave radiation to the snowpack under clear skies, L_dc, is calculated as
 

Ldc = e s Ta4

(9)

where e is effective emissivity, T is average daily air temperature (K), and s is the Stefan-Boltzmann constant. The effective emissivity of atmosphere and canopy overhead, e, is a function of atmospheric water vapor pressure, e_a, atmospheric transmissivity, t, and canopy cover, F. We used Sellers' (1965) coefficients for Brunt's (1932) equation expressing atmospheric emissivity, e_a, as a function of water vapor pressure at dew point temperature. This estimate was adjusted for cloudiness using the equation of Sugita and Brutsaert (1993),
 

C = 1.02 t-0.0227

(10)

where C is a cloudiness scaler and t is atmospheric transmissivity. We used Bristow and Campbell's (1984) procedure to estimate atmospheric transmissivity, as described previously.

The effective emissivity, e, of the atmosphere and canopy overhead is then calculated as
 

e = (1 - F) (0.605 + 0.048 [ea0.5]) C + F

(11)

If the air temperature is below 0 C, then net longwave radiation to the snowpack can be estimated as
 

Ld = (e- 1) s Ta4

(12)

assuming the emissivity of snow as 1.0. When air temperature is above 0 C, the temperature of the snowpack is assumed to be 0 C and net longwave radiation is estimated as
 

Ld = e s Ta4 - s(273)4

(13)

Heat advected to the snowpack by rainfall, R, is calculated based on the assumption that the snowpack is isothermal at the melting point
 

R = rwcw r (Tr - Tm)

(14)

where r_w is the density of water, c_w is the heat capacity of water, r is the amount of rainfall, T_r is the temperature of the rainfall (assumed equal to air temperature), and T_m is the melting point temperature.

All of the energy transferred to the snowpack is assumed to go towards converting snow to liquid (implying an isothermal snowpack) so that 334 kJ are required to melt 1 kg snow.
 

Penman-Monteith equation

Evapotranspiration

Transpiration is calculated using the Penman-Monteith method following a methodology described by Running and Coughlan (1988). The transpiration rate derived from the Penman-Monteith equation m³ H₂0 s day^-1 LAI^-1 is multiplied by daylength (s day^-1) and LAI to obtain canopy transpiration m³ H₂0 s day^-1 for a grid cell.

In BGC the aerodynamic resistance is a constant and assumes a well-ventilated conifer type canopy, an assumption that seems valid for pine and fir forests. However, the model does not distinguish between vegetation types (e.g., grasslands or chaparral) and so modification of the subroutines is necessary to improve characterization of other types of ecosystems. The aerodynamic resistance, r_a = s m^-1, is a reciprocal function of leaf area index (Running and Coughlan, 1988) as implemented in FOREST-BGC
 

ra = 5.0 / (lai / 2)

(16)

which assumes that the canopy is well ventilated with a conifer needle. morphology. This assumption avoids the need for windspeed data.

Table 4 provides a description of the data inputs to the aboveground energy balance subunit of the model and the units.
 

Plant Water Uptake

W_v and W_g are the quantities of water in the vadose and groundwater zones, respectively, T_v and T_g are subject to the constraint that they do no exceed W_v and W_g, respectively.
 

Canopy Conductance

(1) All leaf area receiving a defined minimum threshold radiation, as determined from the canopy light attenuation following a Beer's law extinction equation, will have conductance to water vapor. The Beer's Law canopy light extinction relationship for light environments at depth within the canopy is daily average short-wave radiation,
 

Qi = (Qo(1-ekLAI))/(ekLAI)

(19)

where

Q_i = light transmittance, MJ m^-2 day^-1
Q_o = light intensity above canopy
k = Beer's Law extinction coefficient, assumed to be 0.5
LAI = Leaf area index

In canopies that have light transmission exceeding the radiation threshold, G_c is equal to that defined, in relations (3) and (4) below. If the minimum daily radiation is less than the threshold, then canopy conductance is restricted to cuticular conductance which is defined as 0.00005 m s^-1.

(2) Minimum mean night air temperatures below freezing restrict canopy conductances to the cuticular conductance. Therefore, minimium night temperatures are used to reduce G_c. When below OOC, conductance is reduced by 0.0002 m s^{-1 °} C^-1.

(3) Decreasing soil water supply decreases transpiration through the influence of leaf water potential on canopy conductance. Daily maximum leaf water potential, Y (MPa), is a reciprocal function of soil water supply, where
 

(20)

where

d is soil depth and
f is soil porosity.

Canopy conductance, G_c = m s^-1, is then calculated as
 

Gc = Gc,max - D cc(Y - Y min)

(21)

where

G_c,max is maximum canopy conductance, m s^-1
D _cc is the slope of the canopy conductance versus leaf water potential curve, m s^-1 MPa^-1 (see below),
and
Y _min is the minimum leaf water potential that induces stomatal closure (MPa)
D _cc is calculated from
 

D cc = Gc,max /(Y min - Y min,spring)

(22)

where

Y _min,spring is the spring minimum leaf water potential (Mpa).

(4) Canopy conductance to water vapor is then computed as a function of first leaf water poten . tial and modified by the absolute humidity deficit. Conductance declines under conditions of high vapor pressure deficit, where
 

Gc = Gc,max - slope (Gc,max/Y 1 stomatal closure - Y 1 spring minimum)

(23)

and

G_{c VPD restriction, m s}^-1 = G_{c max} slope (G_{c max, m s}^-1 )/ /Y _{1 stomatal closure, MPa} - Y _{1 spring minimum, MPa}) absolute humidity deficit (m s^-1)

Canopy Interception

Evaporation from Soil

Vegetation Data Sets

Leaf Area Index

Estimation of Leaf Area Index

We also used another method to estimate LAI and % bare soil using a spectral mixing analysis. This technique assumes that the pixel spectrum is composed of the aerial proportional contribution of each of the objects in the pixel. Materials can be distinguished only if they are spectrally distinct. If they have distinctly different spectra then their individual contributions to the mixed pixel can be estimated. This was done for the TM scene and three classes were distinguished: bare soil/litter, green foliage, and "shade/shadow," which captures the surface albedo variation. The contribution of topography to this spectral component was removed by calculating the variation in solar intensity as a function of topography. The residual shade represented a measure of the surface roughness and is related to the spacing and density of the canopy.

Sites lacking vegetation cover (bare soil) were compared to mapped estimates and % cover on the plot datasets.

The spectral mixing equation is
 

Rij = S (Fij Rsoil,j + Fij Rveg,j + Fij Rshade,j) + error

(25)

where

R_ij is the pixel radiance (or calibrated surface reflectance) for each band
F_ij is the fraction of each surface component in the pixel in each band
R_soil,j is the radiance (or reflectance) of soil in each band j
R_veg,j is the radiance (or reflectance) of green foliage in each band j
R_shade,j is the radiance (or reflectance) of shade/shadows in each band j
error is the difference between the R_ij and the sum of the estimated spectral fractions in each pixel

We had no direct measure of LAI or an allometric relationship e.g., using DBH that had been developed for the site and species within our watershed. Comparisons between estimated LAI and canopy greenness were made to % cover from the Forest Service Ecounit site data. A constant scaling factor was applied to convert greenness to LAI using a linear interpolation and scaling between a maximum LAI of 1-12, and is equivalent to maximum NDVI. This comparison provides a estimate of all-sided LAI. Table 5-4 shows that the coefficients of deten-nination (r²) is significant at a probability p>0.0001. The FOREST BGC model uses an empirical relation of 2.2 to divide the all-sided LAI to get pr ecte LAI. For reference, planar leaves in a canopy of spherically distributed leaf angles would have a theoretical value of 2.0.

Table 5-1 shows the results of a comparison between % total cover, % conifer cover, % shrub cover and total % vegetation cover from the Ecounit data set and NDVI and LAI (Table 5-2). None of these relationships had high coefficients of determination (r²) even though some had significant probabilities (P) due to the large sample size (n). This relationship was also tested for topographic slope since some differences in illumination may not have been non-nalized by the NDVI. This relationship also had a low r². We also compared total herbaceous cover, topographic aspect, and hardwood cover against the % shade derived from the spectral mixture analysis after subtracting topographic shade from the analysis (Table 5-3). This relationship also showed a poor fit against the field data sets. The same comparison was performed for the FAI data set with similar results. There is a relatively poor fit between % cover (by forest type) or DBH (trunk diameter at breast height) and greenness as seen by the satellite. This poor fit is probably because of a combination of effors and the fact that as tree cover decreases within a plot, cover by the shrub and herbaceous layer increases. Thus, a better comparison to the satellite data would be total green foliage cover of all vegetation types. Since the model is only parameterized for a single vegetation layer, the greenness probably offers a closer estimate to total plot LAI.
 

Land Cover Type

Soil Properties

Adding more soil layers would be physically more realistic as different species mine different depths of the soil profile or change depth as phenology progresses. The problem is in assigning a root depth (or active root depth) as a function of time, soil type, and vegetation type. There is very little data in the published literature to support a rigorous structuring of rooting depths and root biomass (or other parameters related to water uptake rates, e.g., root length density).

In this study, all of the soil physical and chemical data are from "Soil Survey Eldorado National Forest California", "Soil Survey of Eldorado Area, California", and "Laboratory Data and Descriptions for Some Typical Pedons of CA soil." The soil survey of the western Camp Creek region was done by USDA SCS and the eastern region of Camp Creek region was done by USDA forest service. Both soil surveys are in Order 3 scale. The data used for generating the data file, soilinfo.dat, are taken from Table 5, Table 6, and Table 10 of Soil Survey of Eldorado Area, California (USDA, 1974), and Laboratory Data and Descriptions for Some Typical Pedons of CA soil (SCS, 1973).

The soil attribute data are linked by the soil map units between the soil database and the digital soil survey map. The attribute data for the east side of the study area were directly digitized and input from the hardcopy of the USFS Eldorado National Forest soil survey report. The attribute data for the west side of the study area was abstracted from the STATSCO database.

The attribute data listed in the reports have a range of values (min. and max.) for each soil depth layer. The mean value for each layer was calculated based on the arithmetic average except for permeability which was calculated using the log mean. The model currently treats the soil profile as one layer but could be modified to run at more layers if data is available. The parameters for the soil layer were calculated using a weighted depth. The soil data in the data base includes:

1). Soil type (description)
2). Depth to bedrock (cm)
3). Hydraulic conductivity (cm s^-1)
4). Porosity (dimensionless)
5). Field capacity (dimensionless, cm cm^-1)
6). Effective suction at wetting front (kPa)
7). Available water capacity (dimensionless)

Soil type is a value in the data base linked to the soil map unit symbols (presented by both soil survey agencies). The original soil map units used in the soil survey map and the soil type used in data file can be seen in the data file and are shown in Tables 4 and 5 of both soil survey reports. Depth to bedrock is the sum of the each depth in soil profile and is found in Tables 5 or 6 of the soil survey report.

The effective hydraulic conductivity can be directly estimated from the soil database. However, the saturated hydraulic conductivity which is given in this data base has a wide range; it varies from 0.2 to 6.0 in./hr. for the same soil. Without additional information, we may simply use the average value (or some other algorithm, e.g., log mean, geometric, mean, or harmonic mean, etc.). Saturated hydraulic conductivity was estimated from the original multiple soil layer data set based on the han-nonic mean which is generally used in hydrologic studies for integrating soil moisture. The original data for the west study area can be found in Table 6 of USDA Forest soil survey Soil Survey Eldorado National Forest California.

The porosity of the soil medium can be calculated from the measured bulk density. If the cationexchange capacity of the clay (which is an indicator of the shrink-swell capacity of the clay), and the sand, clay, and organic matter percentage are available, then the bulk density at the water content for 33 kPa tension can be estimated from the following:
 

BD=1.51+0.0025× S-0.0013× S× OM-0.0006× C× OM-0.0048× C× CEC

(26)

in which BD is the bulk density of < 2 mm material (gm/cm³), S is percent sand, C is percent clay, OM is percent organic matter, and CEC is the ratio of cation-exchange capacity (CEC) of clay to the percent clay (CEC ranges from 0.1 to 0.9).

In this study, due to lack of the CEC data, we used sand, clay, and OM to estimation the mineral bulk density following the procedure and data chart developed by Rawls (1985).
 

(27)

Average organic matter bulk density equals 0.224 gm/cm³. Bulk density is converted to porosity accordingly:
 

(28)

where PD is the particle density (2.65 gm/cm³). Specific yield, used in the groundwater model, as assumed to be equal to porosity. The soil particle size distribution data and the chemical data for the east study area are listed in Table 6 and Table 5 of Soil Survey Eldorado National Forest California. Field capacity data are from Table 10 of "Soil Survey Eldorado National Forest California" and "Laboratory Data and Descriptions for Some Typical Pedons of CA soil".

The validity of the soil polygons for use in a hydrologic model of this scale was tested using data from the soil maps which were compared to FS plot data from the Ecounit Inventory. Because soil depth and % clay were the primary infiltration and water holding drivers they were compared by regression analysis. Figure 8a. shows the relationship between the estimated soil polygon % clay content of the soil (weighted average with depth) and the measured % clay content from the Ecounit study. The soil data base did not have a continuous range of % clay contents as seen in the figure. Several Ecounit soil samples were collected for each polygon type and these measured % clay values usually ranged over 20-30%. Thus, even though a significant regression with a positive slope was obtained, considerable scatter is found around the regression line. Not all soils had depth values associated with them but for those that did, we found a significant positive regression between measured depth and polygon value (Figure 8b). Again, substantial scatter is apparent. Nonetheless, these data show that polygon soil properties are within the range of measured values obtained at random point locations. It would be possible to improve the input polygon data by determining a probability density distribution for these soil properties using the field data if this was determined to produce significant error in the infiltration patterns across the landscape.

Effective suction at the wetting front was estimated using the equation given by Rawls and Brakensiek (1982, 1983) by relating the Green-Ampt wetting front suction parameter to soil physical properties in the following equation:
 

(29)

Available water capacity was estimated based on a weighted mean (by the thickness of each layer). The original data for the east side of the study area are listed in Table 6 of Soil Survey Eldorado National Forest California.
 

Digital Terrain Map.

Because topographic gradients drive, both the flow of water and the direction, the model is particularly sensitive to the accuracy of the topographic maps. The resolution of the digital USGS maps are about 30 m in the horizontal plane and 7 m in the vertical plane. We compared the elevation estimated from the USGS topographic data and the derived slope and aspect data to data from the Forest Service Ecounit and FIA datasets. Results from both comparisons yielded similar results with the same conclusions. Data from the Ecounit study are shown in Figure 9. It shows that the elevation predicted from the DEM is very close to the measured elevation. Aspect is also relatively well predicted from field measurements, while the worst case is found for estimating slope. The greater scatter in the predicted slope apparently results from both greater problems in mapping it accurately at the field scale and the differences in scales measured. When DEM data are aggregated to 5 x 5, 9 x 9, and 15 x 15 pixels, the relationship between measured and predicted values decreases. Thus, slopes are under-estimated which affects the accuracy in predicting runoff for all hydrologic models.

The DEM for the Camp Creek watershed is shown in Figure 10. The smaller Snow Creek watershed is shown in Figure 11 and the Pebble Creek watershed in Figure 12. Snow Creek and Pebble Creek are located within the watershed shown in Figure 10. The location of Snow Creek within the Camp Creek watershed can be seen in Figure 1.
 

Synthetic Terrain Map.

The model was tested by several sets of generated surface flow parameters:

- Three rainfall types : point, polygon or constant over the area
- Three constant slopes : 10, 30, or 60 %
- Rotated in two kinds of grid orientations : normal to the N or to the NE direction
- Different shapes of landscape : steep plane, divergent, and convergent steep surfaces

- Constant slope
- Composed of four valleys having four different slopes
- Octagonal DEM; created to test anisotropy
- A cone
- A valley with converging and/or diverging surfaces

Surface Flow

In the first case the soil is saturated, in the second case the precipitation rate exceeds the infiltration rate, and in the third case exfiltration from groundwater contributes to the surface flow.

Knowing precipitation, infiltration, and discharges from the adjacent cells, a change in surface water storage in the cell of interest is:
 

(33)

in which P is precipitation rate, I is infiltration rate, D x and D y are the cell dimensions in the x and y directions, respectively, and Q_s is the flow from the four adjacent cells during the time step D t. Discharge is the product of velocity and area. Flow area is D x or D y multiplied by flow depth h_s. If inflow exceeds outflow, the h_s increases or if more water leaves than enters, it falls.

Water moves across adjacent cell boundaries based on physical processes and the topographic conditions. Figure 5 shows how the surface flow is calculated based on a two dimensional representation of water flux. The slope of the water surface is assumed to equal the slope of the ground surface, which is usually acceptable for gradually varied steady flow. Overland flows accumulate in depressions and in small streams which eventually merge into rivers, lakes, or ponds. In this phase, the transport of the water can be considered as free surface, gradually varied, sheet flow. The energy source for surface flow is gravity, which is consumed by friction, so the overland flow is considered as gradually varied sheet flow. Manning's equation states that the overland flow velocity -is a function of surface geomorphologic and surface water ponding conditions. For the x direction,
in which u is velocity m/s, n is Manning's friction coefficient L^-1/3 T^-1, R is hydraulic radius m (flow section area A m² over the wetted perimeter WP m), and S is bed slope m/m. For sheet flow, area is just flow depth h_s multiplied by unit width and, thus, hydraulic radius simplifies to flow depth.
 
 

(34)

in which u is velocity m/s, n is Manning's friction coefficient L^-1/3 T^-1, R is hydraulic radius m (flow section area A m² over the wetted perimeter WP m), and S is bed slope m/m.  For sheet flow, area is just flow depth h_s, multiplied by unit width and, thus, hydraulic radius simplifies to flow depth.
 

	(35)
Qx = uAx	(36)

and similarly for the y direction,
 
 

	(37)
Qy = vAy.	(38)

There are two faces for the x direction and two for the y direction for a total four, as given in the mass balance equation Q_s.

Appropriate use of Manning's equation assumes particular conditions. Velocity is constant with depth; flow is uniform, that is depth, wetted section, velocity and discharge are constant along the channel; and the slope of the water surface is equal to the linear loss of head along the channel. For gradually varied flow, the equation applies at a location but the velocity changes as slope, roughness and depth varying spatially over the landscape. A constant Manning's n (0.1) was used because litter depth did not vary with vegetation types in this data base.

Field measurements (Starosolszky, 1987) showed that the Manning's equation yields estimates which agree with measured values for water at typical environmental temperatures (about 20° C) and flows at moderate depth (< 10 meters). Using this method for estimating the surface flow simplifies the model and minimizes parameter input into variables more readily derivable in GIS applications.
 

Subsurface Flow

Lateral subsurface flow through the vertical edges or faces of each element is controlled by the thickness of the saturated layer and the slope or gradient of table elevation. There is no leakage through the bedrock. Since the slope or change in the depth of the water table is small between adjacent cells and the unconfined flow field is shallow in our system in the absence of sources and sink, the Dupuit assumption is satisfied. The Dupuits' assumption means that the flow moves in an horizontal plane (the isopotential lines are vertical). Darcy's law is limited to:
 

steady flow, ¶ h/¶ t = 0 and laminar flow.

(40)

Otherwise energy is wasted by viscous friction. But usually the following limits are acceptable: R_e represents the Reynolds' number, or the ratio of inertial forces over viscosity forces.
 

if Re <1,	flow is laminar and Darcy's law applicable.
if 1 < Re <60,  qn = -K D H/L	Darcy's law is applicable with the correction proposed by Dupuits:
if 60 < Re	Darcy's law is not applicable.

In subsurface flows the velocity is usually very slow and Darcy's law can be used without considering the Reynolds' number. There is a problem with Darcy's law application for verv low and very high values of the hydraulic conductivity K. Conductivity has no meaning in the case of a cracked rock, for example. In this case study, the conductivity at saturation ranges between 7.10^-1 and 2.8 10^-4 ms^-1.

To simplify the problem, instead of solving Poisson's equation for steady state flow, Darcy's equation was solved directly for subsurface flow. The discharge Q per unit width for an aquifer thickness b in the x and y directions, respectively, is:
 

(41)

in which k is saturated hydraulic conductivity; hydraulic gradient is the change in total head with distance
in the x and y direction and , respectively. Total head in each cell is the elevation head defined by the base of the active soil layer plus the water table thickness b and the surface ponding depth h_s, if the water table is at or above the ground surface. The datum for the elevation head is the same as that used for the interpolated DEM data. Spatial distance along the head drop between two cells is the grid spacing.
 

Vadose zone

The vadose zone mass balance equation is coupled to the groundwater equation because h_v, depends on h_s, (groundwater elevation). To decouple the mass balance equations, the thickness of the vadose zone h_v is assumed constant during a time step and unadjusted water content q * is calculated. At the end of the time step, depth of the vadose zone and water content are adjusted to reflect the change in groundwater elevation ().
 

	(43)
	(44)

The water content of the incremental layer is field capacity.
 
 

Infiltration

In the case of no surface water ponding, infiltration rate is limited by the precipitation rate (model input). With ponding, infiltration rate is controlled by surface or subsurface conditions. Surface conditions constrain infiltration according to the potential infiltration rate I_ga estimated using the Green and Ampt equation (1911) or according to the depth of water ponding on surface which could infiltrate during the time step. Subsurface conditions might constrain infiltration according to the amount of unfilled storage in the vadose zone. If the vadose zone is filled, the groundwater table is the free water surface and the groundwater system can contribute to (exfiltration) or gain water from the surface flow.

Almost all subsurface flow models are based on Darcy's Law, and most approximate models give similar results. Vertessy et al. (1993) used a ID Richard's equation to solve for unsaturated flow and a 2D Richard's equation was used for saturated flow over a 0.32 km² forest watershed. However, this detailed physical approach has not been applied to a large watershed.

The Green and Ampt (1911) model is most frequently used for larger watersheds because model assumptions simplify parameter estimation. In the Green-Ampt equation, the initial infiltration rate, the saturated infiltration rate, and soil physical properties are considered. The original Green-Ampt infiltration model was developed for ponded infiltration into a deep homogeneous soil with a uniform initial water content. In applying the Green-Ampt approach, two assumptions are made. First, there is a distinct and precisely definable wetting front. The matric potential at this wetting front remains effectively constant, regardless of time and position (i.e., dependent only on soil physical properties). The second, is that behind the wetting front the soil is uniformly wet and of constant conductivity. By neglecting the surface ponding depth the Green-Ampt rate equation has the form:
 

(45)

in which I_ga is potential infiltration rate m/s, k is effective hydraulic conductivity, m/s, H_f is effective matric potential at the wetting front m, M_d is the difference between saturation and initial water content (f -q _i) m/m, q _i is water content at the beginning of the precipitation event m/m, f is soil porosity m/m, F is cumulative infiltration m. Notice that the infiltration rate is dependent on the water content at the beginning of the precipitation event which is determined by the vadose zone mass balance equation. In forward finite difference form,
 

	(46)
	(47)
	(48)

in which

a = 1

Mein and Larson (1973) adopted the Green-Ampt model for infiltration during rainfall. Just prior to surface ponding the infiltration rate I equals the rainfall rate R and the cumulative infiltration at the time to surface ponding is the same as the cumulative rainfall. Over a short time period (e.g., 1 second), we assume the rainfall is homogeneous and steady in each cell. The available water for infiltration is deten-nined for each cell by the contributions of both rainfall and lateral surface flow. Ponding occurs when available water exceeds the sum of the rainfall and lateral surface flow. The surface water not only comes from rainfall, but also from adjacent cells via surface flow, so the cumulative infiltration is determined by both infiltration rate and the water supply from the surface.

The soil column becomes saturated when cumulative infiltration is greater than or equal to f -q _i times the soil depth. Next, the infiltration rate is limited by saturated hydraulic conductivity. When the subsurface is saturated, cxfiitration may occur due to the subsurface flow. The surface-subsurface flow system is dominated by changes in the hydraulic head and changes in subsurface storage.
 

Recharge

R=0 if I ¹ 0
R=k if I=0 and q _fc £ q £ f

Conversely, if infiltration stops before the advancing front meets the water table there is recharge. Water in the vadose zone, retained between field capacity and saturation and behind the wetting front, drains or recharges the groundwater at the rate of the effective hydraulic conductivity. Evapotranspiration is removed from the vadose zone if the water table is below the soil surface and it is removed from the groundwater flow compartment when the watertable rises above the soil surface.
 

Evapotranspiration

Boundary and Initial Conditions

We used the assumption of Famiglietti and Wood (1994b) for initializing the hydrologic model by assuming that the local soil profile state is at gravitational equilibrium. Our model only used a one layer soil and we assumed that at the beginning of the simulation that there was no surface water ponding and that the total head was equivalent to the elevation head (i.e., there was no water table) except where free surface water existed. For the lakes, river channel, and streams, the total head is the sum of elevation head and the hydraulic head. This definition explicitly defines the surface and subsurface flows which are zero initially.
 

SNEP Scenario Testing and Simulation Results

The most interesting observation is the large spatial pattern in ET within this watershed (as observed by the gray scale where black = lowest values and white equals highest. The maximum annual ET and weighted % vegetation cover is shown in Table 6. Transpiration was found to increase by 10.98% over the five decades while cover increased by 11.55%. To perform this analysis we disaggregated the 74 polygon cover values to 3939 pixels while the HRU model aggregated these patterns into three HRU polygons. The level of resolution is apparent in Figure 16 which shows the LAI derived from the NDVI analysis for this area. Several observations can be made, illustrating the value of this model. First, there is more continuous variation in LAI in the satellite data. Second, although the satellite patterns are generally observed in the polygon data, there is a difference in both location and size of the polygons relative to the satellite imagery. While the exact LAI values of the NDVI may be questioned, the spatial pattern represents more correctly the actual vegetation patterns than the polygon representation. Lastly, the initial LAI values derived from the satellite and shown in Figure 16 are higher than those in Figure 15, which were based on an allometric relationship. The difference in starting LAI at decade 0 is substantially greater (24% higher, assuming a maximum allsided LAI of 12) than that for the 50 year change used in this scenario (11.55% difference over the 50 years). This result indicates the sensitivity of these models to input variables and their accuracy.

Results from the surface-subsurface model include: spatial patterns of cumulative precipitation, infiltration rate, surface flow velocity, subsurface flow velocity, water table elevation, surface water ponding depth, vadose zone depth. The rainfall event started at time zero and continued for seven hours. Precipitation increases with elevation, following the DEM (Figure 11).

The effect of spatially and temporally varying rainfall rates as well as spatially varying soil types and initial soil moisture contents, produce the temporal and spatial patterns observed in the instantaneous infiltration rate (not shown). The infiltration rate and cumulative infiltration are important factors in understanding pollution management. For example, essential soil nutrients are largely lost through vertical leaching. Thus, spatial maps of infiltration are of significant importance in understanding biogeochemical cycling in boreal ecosystems beyond the immediate role in the partitioning of the hydrologic budget.

The presence of surface water ponding and its spatial and temporal pattern during and after rainfall events is an important parameter for land use, flood control planning, and for hydraulic engineering. Figure 17 shows spatial patterns in surface ponding at 5 hours. From this image, one observes the strong dependence of surface pending on terrain features. Surface poinding depth is 6 mm at the top of the hills and is 4 mm in the trough after 5 hours. Although not shown surface ponding depth in sites having steep slopes and high elevation has decreased, but the depth in the trough continues to increase. Later, most surface ponding is restricted to stream channels. The spatial and temporal patterns in surface ponding depth demonstrate the effect of both terrain features and precipitation patterns on surface water storage.

The surface flow velocity is important in the study of soil erosion, sediment yield, sediment transport, hydraulic engineering designs, soil nutrition, and redistribution of water resources. Figure 18 shows the spatial patterns in surface flow velocities at 7.4 hours. In the mountain regions, where the terrain is steep and the rainfall rate is high, the surface flow velocities are higher than in the trough. The surface flow velocity is equal to or greater than I m/s on steep slopes and is only a fraction of this in the low gradient trough areas.

Subsurface flow velocity is important for contamination studies as well as characterizing water resources distribution. The water table surface above DEM datum (Figure 19) demonstrates the variations in subsurface water resources. The water table reaches the soil surface at the time when the soil column becomes completely saturated. Drainage of water from the surface to the depth of the active layer is controlled by subsurface flow and infiltration. When the water table rises, the water table comes closer'to the ground surface. At 7.4 hours the water surface is higher than in the trough and this gradient drives the groundwater to the trough. Because the gradient is steeper at higher elevation, the subsurface velocity is higher (Figure 20). Although the soil type varies (Figure 21), the conductivity is the same except for the south west comer. Soil in this comer is deeper and has a greater hydraulic conductivity than in remainder of the watershed. The water table drops sharply as shown by the spike in the thickness of the vadose zone (Figure 22, but not obvious from Figure 19 due to scale of the z axis).

 

REFERENCES

Beven, K. J. and M. J. Kirkby. 1979. A physically based, variable contributing area model of basin hydrology. Hydrologic Science Bulletin. 24:43-69.

Bristow, K. L. and G. S. Campbell. 1984. On the relationship between incoming solar radiation and daily maximum and minimum temperature. Agricultural and Forest Meteorology. 31:159-166.

Brunt, D. 1932. Notes on radiation in the atmosphere. 1. Quarterly Journal of the Royal Meteorological Society. 58:389-420.

Courant, R., K. 0. Friedrichs, and H. Lewy. 1928. Uber die partiellen differenzengleichngen der mathematishen physik. Mathematishe Annalen. 100:32-74.

Dickinson, R. E., A. Henderson-Sellers, C. Rosenzweig, and P. J. Sellers. 1991. Evapotranspiration models with canopy resistance for use in climate i-nodels. Annual Review: Agricultural Forestry Meteorology. 54:373-388.

Eagleson, P. S. 1978a. Climate, soil, and vegetation, 3. A simplified model of soil moisture movement in the liquid phase. Water Resources Research. 14:722-730.

Eagleson, P. S. 1978b. Climate, soil, and vegetation, 4. The expected value of annual evapotranspiration. Water Resources Research. 14:731-739.

Famiglietti, J. S. and E. F. Wood. 1991. Evapotranspiration and runoff from large land areas: land surface hydrology for atmospheric general circulation models. Survey in Geophysics. 12:179-204.

Famiglietti, J. S. and E. F. Wood. 1994b. Application of multiscale water and energy balance models on a tallgrass prairie Water Resources Research. 30:3079-3093.

Gao, X., S. Sorooshian, and D. C. Goodrich. 1993. Linkage of a GIS to a distributed rainfall-runoff model. In Environmental Modeling with GIS. Edited by M.F. Goodchild, B. 0. Parks, L. T. Steyaert, 182-187. New York: Oxford Press.

Glassy, J. M. and S. W. Running. 1994. Validating diurnal climatology logic of the MT-CLIM model across a climatic gradient in Oregon. Ecological Applications, 4:248-257.

Green, W. H. and G. A. Ampt. 191 1. Studies on Soil Physics: 1. Flow of air and water through soil. Journal of Agricultural Science. 4:1-24.

Hungerford, R. D., R. R. Nemani, S. W. Running, and J. C. Coughlan. 1989. MTCLIM: A mountain microclimate simulation model. In Research Paper INT-414. Ogden, Utah: USDA For. Ser. Intermountain Research Station.

Jordan, J.P. 1992. Identification des processus de gration des crues -Application au bassin de la Haute Mentue" - Ph.D dissertation. Ecole Polytechnique rale de Lausanne, Switzerland.

Kattelmann, R. 1990. Effects of forest cover on a snowpack in the Sierra Nevada. In Watershed Planning and Analysis in Action. edited by R. E. Riggins, E. B. Jones, R. Singh and P. A. Rechard. 276-284. American Society of Civil Engineers. New York: New York.

Kattelmann, R.C., N. H. Berg, and M. K. Pack. 1985. Estimating regional snow water equivalent with a simple simulation model. Water Resources Bulletin 21(2):273-280.

Laramie, R.L. and J. C. Schaake, Jr. 1972. Simulation of the continuous snowmelt process. In Report 143. Cambridge, Massachusetts: Ralph M. Parsons Laboratory, Massachusetts Institute of Technology.

Leavesley, G. H., R. W. Litchy, B. M. Troutman, and L. G. Sandon. 1983. Precipitation-runoff modeling system: User's manual. In U.S. G. S. Water Resources Investigations Report 83-4328. p. 209.

Maidment D.R. 1993. Handbook of hydrology. New York: McGraw-Hill.

Mein, R. G., and C. L. Larson. 1 973. Modeling infiltration during a steady rain, Water Resources Research. 9:384-394.

Monteith, J. L. 1963. Gas exchange in plant communities. Environmental control ofplant growth. edited by L. T. Evans. 95-112. New York: Academic Press.

Penman, H. L. 1948. Natural evaporation from open water, bare soil, and grass. In Proceedings of the Royal Society of agriculture. 193: 120-145.

Rawls, W. J. and D. L. Brakensiek. 1983. A procedure to predict Green and Ampt infiltration parameters, Advances in Infiltration. 102-112. American Society of Agricultural Engineers.

Rawls, W. J., D. L. Brakensiek, and K. E. Saxton. 1982, Estimation of soil water properties, Transactions of the American Society of Agricultural Engineering. 25:1316-1320, 1328.

Rawls, W. J. and D. L. Brakensiek. 1985. Prediction of soil water properties for hydrologic modeling. In: Watershed Management in the Eighties American Society of Chemical Engineers. 293-299.

Running, S. W. and J. C. Coughlan. 1988. A general model of forest ecosystem processes for regional applications. 1. Hydrologic balance, canopy gas exchange and primary production processes. Ecological Modeling. 42:125-154.

Running, S. W., R. R. Nemani, and R. D. Hungerford. 1987. Extrapolation of synoptic meteorological data in mountainous terrain, and its use for simulating forest evapotranspiration and photosynthesis. Canandian Journal of Forestry Research. 17:472-483.

Soil Conservation Service. 1985. Soil Survey of the Eldorado National Forest California. Pacific Southwest Region: U.S.D.A. Forest Service.

Soil Conservation Service. 1973. Soil Survey laboratory data and descriptions for some soils of California. Soil Survey Investigations Report No. 24. Washington, D.C.: U.S.D.A. Soil Conservation Service.

Sellers, P. J., Y. Mintz, Y. C. Sud, and A. Dalcher. 1986. A simple biosphere model (SiB) for use within general circulation models. Journal of Atmospheric Science. 43:505-531.

Sellers, W.D. 1965. Physical Climatology. Chicago, Illinois: University of Chicago Press.

Sinclair, TA.R. 1994. A diagnostic model for estimating orographic precipitation. Journal of Applied Meteorology. 33:1163-1175.

Smith, J.L. 1982. The historical climatic regime and the projected impact of weather.

Starosolszky, 0. eds. 1987. Applied Surface Hydrology. Littleton, Colorado: Water Resources Publications. 40-120.

Sugita, M. and W. Brutsaert. 1993. Cloud effect in the estimation of instantaneous downward longwave radiation. Water Resources Research. 29(3):599-605.

USDA Soil Conservation Service. 1974. Soil Survey of Eldorado area, California. Washington, D. C.: U.S.D.A. Soil Conservation Service and Forest Service.

Vertessy R. A., T. J. Hatton, P. J. O'Shaughnessy, and M. D. A. Jayasuriya. 1993. Predicting water yield from a mountain ash forest catchment using a terrain analysis based catchment model. Journal of Hydrology. 150:665-700.

Wigmosta M. S., L. W. Vail, and D. P. Lettenmaier. 1994. A distributed hydrology-vegetation model for complex terrain. Water Resources Research. 30:1665-1679.

Xiao, Q. F., S. L. Ustin, and W. W. Wallender. 1996. A spatial and temporal continuous surface-subsurface hydrologic model. Journal Geophysical Research-Atmospheres. 101(D23):29565-29584.

Xiao, Q. F. 1994. A spatially and temporally continuous surface-subsurface hydrologic model. Master's thesis, Department of Land, Air, and Water Resources, University of California, Davis.