Hydrologic Science Section Department of Land, Air, and Water Resources University of California, Davis, CA 95616   Revised Manuscript submitted to: Journal of Geophysical Research-Atmospheres July 1995

Abstract

INTRODUCTION

For many applications, detailed spatial and temporal information on hydrologic processes is needed.  Satellite remote sensing data can provide the detailed spatial information on which to drive spatially varying hydrologic models.  Further, satellite sensors can provide not only synoptic coverage at high ground resolution but also provide repeated sampling over time.  Satellite sensors collect data at intervals from days to months, in various spectral bands at wavelength regions that may directly measure some environmental conditions and provide additional sources of assessment information.  Geographic Information System (GIS) software has powerful spatial data analysis capabilities and provides a means to integrate traditional hydrologic models with spatial data inputs, e.g., from satellite data sources.

In recent years, remote sensing and GIS techniques have become more widely used in hydrological and ecological studies (Maidment, 1993; Moore et al., 1993).  GIS applications have ranged from water resources management to ground water transport (Nachtnebel et al., 1993).  For example, Kaden (1993) used GIS to examine nitrate transport in soils and groundwater contamination and Vieux (1991) modeled non-point source water quality and quantity.  A ground water model was developed by Batelaan et al. (1993).  In most cases, surface and subsurface transport processes have not been linked or models have not been physically based, although recently new models have addressed these issues.  Goodrich (1990) and Gao et al. (1993), developed a fine-scale (15 m) spatially explicit physically based surface-subsurface model for a small  (4.4 ha) semiarid watershed using a digital elevation model (DEM), soil, vegetation, and land use maps.  They used a 1D Richard’s equation to simulate unsaturated flow, a 2D groundwater flow model that assumed horizontal flow, and a 2D kinematics wave equation to simulate overland flow.  Stream flow was routed using the St. Venant equation. It was possible to model this basin implicitly by numerically solving partial differential equations because the area was small and intensively instrumented, and since groundwater flow was deep and isolated from the surface, it could be ignored.  Although this model was physically robust, it is impractical to apply at the scale of typical watersheds because the data needed to drive the model is unavailable.

Use of satellite images to parameterize hydrologic models and provide up-to-date information on surface conditions has also become more common.  Systeme Probatoire d'Observation de la Terre (SPOT, 20 m resolution), Landsat Thematic Mapper (TM, 30 m), and the NOAA Advanced Very High Resolution Radiometer (AVHRR, 1.1 km and 4.4 km) satellite remote sensing data have been widely used for ecological studies.   For example, Band (1993) used AVHRR data to study the impacts of spatial scale and representations on hydrologic processes and Ostendorf and Reynolds (1993) used SPOT data to compare vegetation distribution to runoff patterns they observed in a 2.2 km² Imnavait Creek watershed on the northern foothills of the Brooks Range, Alaska.

A number of recent studies have examined hydrologic processes under extreme climate conditions found in arctic and arid lands.  Ostendorf and Reynolds (1993) developed an elevation-driven surface flow model called T-HYDRO that transported water between grids by assuming pipe flow. A Swedish conceptual hydrologic runoff model, HBV, was developed by Bergström (1976) and used by Hinzman and Kane (1991, 1992) in their extended study of hydrology of the Imnavait Creek basin. This model was run at one-day time steps and generated runoff by summing subbasin contributions.  The subbasins were defined by elevation zones and vegetation type which were found to significantly impact calculated runoff.  The lumped parameters used in the HBV model are very similar to HRU’s (hydrologic response units) used by Ross et al. (1979) in the Finite Element Storm Hydrograph Model (FESHM), Leaversley et al. (1983) in the Precipitation-runoff modeling system (PRMS), or to the soil-vegetation index used by Band (1993) and Band et al. (1993).

More recently, the catchment-scale water balance model of Famiglietti and Wood, (1994a), the coupled hydrological-ecological model of Band et al. (1993), or the distributed hydrology-vegetation model of Wigmosta et al. (1994) were developed to study linked ecological and hydrological processes.  Models such as the Regional HydroEcological Simulation System (RHESSys) (Band, 1993; Band et al., 1993; Nemani et al., 1993), and the Soil-Vegetation-Atmosphere Transfer Scheme (SVATS) by Famiglietti and Wood (1991, 1994a,b, 1995) combine a set of remote sensing and GIS techniques with an integrated spatially distributed hydrologic model.  The RHESSys model is derived from TOPMODEL (Beven and Kirkby, 1979), a quasi-distributed hydrologic model which is linked to the FOREST-BGC model (Running and Coughlan, 1988), an ecosystem process model predicting carbon and water budgets.  The statistically based TOPMODEL (also used by Famiglietti and Wood) differs from earlier soil-vegetation process models by including spatial variations in the main driving processes and uses topographic HRU’s to calculate lateral soil water transport.  Nemani et al. (1993) used AVHRR images to estimate LAI and ran the RHESSys model at daily time steps and two levels of spatial aggregation, based on HRU’s defined by elevation and LAI, over the Seeley-Swan valley (defining 5 or 170 HRU partitions within the 5000 km² region)  in northwestern Montana.

The HydroEcological simulation for Nemani et al. (1993) and Band (1993; Band et al. (1993), was derived from weather station climate data that was modified by local elevation using a weather simulator, MT-CLIM (Running et al., 1987).  Famigletti and Wood (1991, 1994a,b, 1995) and Famigletti et al. (1992) examined spatial variability in water balance and evapotranspiration from the King’s Creek catchment as part of the Konza Prarie, Kansas FIFE experiment.  The model evaluated diurnal dynamics of water and energy fluxes using half-hour time steps and DEM resolved at 30 m.  The SVATS assumed no lateral flow in the unsaturated soil zone and runoff and ET patterns were statistically represented. These distributed or quasi-distributed model studies are primarily aimed at the subbasin scale of hydrologic processes (smaller than 25 km²).

Wigmosta et al. (1994) describe a physically based distributed hydrology-vegetation model they developed for complex terrain and tested on a 2900 km² watershed in northwestern Montana. This model was run at 3 hr time steps with 180 m spatial grid and links spatially explicit dynamic patterns of surface energy budget, runoff, subsurface flow and soil storage.  The model was calibrated by adjusting soil properties to fit a two year precipitation and stream flow data record and tested against two additional years of data.  No overland flow or infiltration were computed in the model.

Vertessy et al. (1993) also developed a distributed parameter terrain model, Topog_yield, which simulates slowly changing hydrologic states over multi-year sequences.  Their model, which was run at daily time steps for 12 years over a 0.32 km² catchment showed only a small cumulative runoff error compared to observed stream flow.  Even though these examples and other hydrologic models have been developed using GIS techniques and remote sensing data, most are better parameterized to represent seasonal and inter-annual processes (e.g., for understanding plant water use) and for studying long-term hydrologic processes then for short-term hydrologic event predictions.

For understanding storm events and for flood forecasting more rapidly changing hydrologic processes need to be considered and analyzed over large watersheds.  Most models are limited by their spatial resolution. Typically, spatially defined hydrologic models either have course resolution (e.g., 1 km) or are well instrumented but small catchments of a few ha. Models applied to larger watersheds are generally empirically calibrated using hydrograph records while physically based models have been largely restricted to small catchments.  In most cases, there is no linkage between the surface flow, soil storage, and subsurface flow processes.  Additionally, many of these models run at daily to monthly time steps and all processes are defined at the same temporal resolution.  These scale problems introduce artifacts as the models and the processes are not operating at the same space/time scales.  In fact, the scaling issue is a major research focus today (e.g, Band, 1993; and Famiglietti and Wood, 1994, 1995).

In this research, a cell-based hydrologic model linking transport between surface and subsurface hydrologic process is built based on continuity theory and mass balance.  All of the parameters needed by this model are directly or indirectly derived from satellite remote sensing data, USGS topographic maps, NRCS soil maps, climatological and hydrologic data, and parameters obtained from the literature.  This model is based on physical relationships, and as such, it should be generally suitable for application in other areas.  The time-variable model is linked to the GIS primarily through the use of database management and visualization tasks.

MODEL DESCRIPTION

The study area is gridded at the same scale in X and Y directions.  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.  Figure 2 illustrates how the grid and the model are abstracted from the natural landscape.  A mass balance is computed for each cell in the grid.  Figure 3 (a) shows a three dimensional representation of the hydrologic processes in the model.  Figure 3 (b) shows a X- Y plane view of the surface and subsurface flow paths.

We applied our model to a large watershed in the Yukon River Basin near Circle, Alaska.  This site was chosen to develop and test the event based model because a number of simplifying assumptions were justified.  First, we modeled a single summer convective storm event (August 14, 1987) and the hydrologic processes after this event.  This time period was selected because precipitation input occurs as rainfall rather than snow and because precipitation occurs in isolated storm events (lasting a day to a few days) which permits the active layer to drain between storms.  Relatively frequent storm events occur in the watershed during July, August and September.  Kane et al. (1989, 1990) describe similar seasonal storms for the north slopes of the Brooks Range.  In midsummer, the snow pack has melted and the active layer above the permafrost has melted to near maximum depth.  The surface layer of organic soil (below the litter) remains relatively wet, so the water volume absorbed by wetting the organic layer was neglected.  In this region, vertical subsurface movement of water is limited by the depth of permafrost.  The permafrost depth is determined by the vegetation cover type, density of the canopy cover, topography, and soil type following the conceptual model of Van Cleve et al. (1983).

A further assumption in this model exercise is that because the humidity is high during a storm event, surface evapotranspiration could be ignored for the duration of the storm and for a few days thereafter. There is little data on evapotranspiration rates in Alaskan forests, although evapotranspiration has been estimated to consume 32% and  34%-66% of the annual precipitation in central Alaska (Anderson, 1970) and arctic tundra watersheds (Kane et al., 1990), respectively.  Patric and Black (1968) computed potential evapotranspiration using Thornthwaite’s equations (1948) and estimated 350-450 mm/yr for interior Alaska sites with precipitation deficits of 290 mm at Fort Yukon (60 km north of our watershed) and 190 mm at Fairbanks (180 km south of our watershed).  If the growing season is assumed to be four months long, the mean daily ET is 1-2 mm.  Dingman (1966) estimates that ET consumes 76% of the summer precipitation near Fairbanks, and Kane et al., (1980, 1990) also estimate summer evapotranspiration of only 1.0 mm/day for bare soil to 2.0 mm/day for tundra on the Alaska north slope.  This level of potential daily ET loss over taiga and tundra sites suggests that losses over our watershed are much smaller than the cumulative precipitation of a single large summer storm.  Thus, errors induced by not accounting for ET during and for a few days immediately after a summer storm are relatively minor.

Finally, a characteristic of the boreal forests of the region is that trees are small and canopies are sparse.  These low stature forests with open canopy structure allowed us to assume that interception flow within the overstory canopies was minimal and could be ignored for a first approximation. The understory, particularly moss, grasses, and sedges can hold substantial amount of water so we assumed that the ground surface was wet and free evaporation could occur but that the total amount was low because of low net radiation.  Because there was a free path for water flow into the soil, interception by the understory was ignored. These assumptions simplified the model for calculating the hydrologic budget over a one-week period.
 

Precipitation

For calculating the rainfall rate, we divided the study region into three elevation zones: flood plain, old river terraces (lowlands), and montane (Xiao, 1994) while assuming that the prevailing winds flow from south to north as determined from the Circle City NOAA weather records.  Based on the empirically derived estimates of Zhang et al. (1987) the rainfall rate was estimated by:
 

	P = cP0[a(E - E0) + b(b - b0)] + P0	(1)
where	b  = a - 90       "       0 </= a< 180
	b = 270 - a      "        180 </= a </= 360

in which P is rainfall rate adjusted for elevation and aspect.  E is elevation (m), E₀ is the weather station elevation (182.3 m for this case study),  b is adjusted site aspect, b₀ is adjusted aspect of the weather station (0 for this case study), a is true aspect,  P₀ is the instantaneous rainfall rate (mm/hr) at the weather station, a is the change of rainfall rate with elevation.  For this case study, the change with elevation was equal to 33 mm / 100 m * (time step * P₀ )/(mean total annual precipitation), and b is the change of rainfall rate with aspect which is 40 mm / 180^o  *(time step * P₀ )/(mean total annual precipitation), and c is a coefficient which is used to remove micro-topographic effects on precipitation rate.  In this study c = 1 for mountain uplands,  c = 0.75 for terraces, and c = 0 for flood plains.
 

Surface Flow

Three types of surface runoff are observed: saturation excess runoff, infiltration excess runoff, and subsurface flow that are described in detail by Famiglietti and Wood (1991).  In the first case, runoff occurs immediately when precipitation falls on saturated soil, in the second case, runoff occurs because precipitation exceeds the infiltration rate and, in the third case, runoff from adjacent areas occurs when subsurface flow exits the hillslope. The model captures all three phenomenon and Manning's equation is used to calculate overland flow velocity as a function of surface geomorphologic and surface water ponding conditions:
 

(2)

in which n is Manning's surface roughness coefficient and R is the hydraulic radius, the ratio of water cross-sectional area to its perimeter.  For sheet flow, R simplifies to the surface water ponding depth, and S is the surface slope.

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.

Surface ponding time was defined from the rainfall inception and the depth of surface ponding and these variables became important parameters for the surface flow sub-model.  At the beginning of the rainfall event, due to the low initial water content of the soil, the theoretical infiltration rate was higher than the rainfall rate.  With time, cumulative rainfall increased and infiltration increased with rainfall.  Finally, at some time step, the rainfall rate exceeded the infiltration rate, and the excess rainfall started ponding, at which point, surface flow started.  Those areas which have local topographic sinks receive the surface flow until the sink fills and flow spills over into adjacent cells.

Figure 3 and Equation (2) show how the surface flow is calculated based on a two dimensional representation of water flux.  We assumed that at each time step, the flow into or out of a cell is only influenced by adjacent cells.  Only eight flow directions are possible, the cardinal directions and their intermediate diagonal directions (e.g., southeast).  Flows to the north, east, south, and west can be directly calculated from Equation 2.  For this case, the cross-section area of the surface water flux is the product of the cell size and the surface ponding depth.  To calculate flows in the diagonal directions requires computing the proportional velocity in the x and y directions.  For water flow in the diagonal directions, the flow velocity in the x-direction is the product of the cos 45º * Vmean, and the flow velocity in the y-direction is the product of the sin 45º * Vmean.  The cross-section area surface flow is the same as in the previous case, since we use equal cell sizes in both directions.

The Manning’s surface roughness coefficient is estimated from land use, vegetation type, and surface cover.  We use the surface water ponding depth as the hydraulic radius, and the surface slope is determined from the DEM.  For the zero slope case,  instead of using only the DEM data to calculate surface slope, the surface water ponding depth is used to calculate the hydrologic slope so that a smoother and more rational surface flow pattern is obtained.
 

Infiltration

The Green and Ampt (1911) model is most frequently used because assumptions simplify parameter estimation.  In the Green-Ampt equation, the initial infiltration rate, the saturated infiltration rate, and the 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 principal 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:
 

(3)

and the integrated form is:
 

(4)

where K is the effective hydraulic conductivity;  S_f  is the effective matrix potential at the wetting front; f is the soil porosity; q_i is the initial water content; F is the cumulative infiltration; and f is the infiltration rate.

Mein and Larson (1973) adopted the Green-Ampt model to infiltration during rainfall.  Just prior to surface ponding the infiltration rate f equals the rainfall rate P and the cumulative infiltration at the time to surface ponding is the same as the cumulative rainfall. In this study, a lateral surface flow velocity from adjacent cells is considered and is presented later.   The available water for infiltration is determined for each cell by the contributions of both rainfall and lateral surface flow.  The relationship among rainfall rate P,  lateral surface flow velocity from adjacent cells, and the infiltration rate f is described as:
 

f = P + net lateral surface flow for t </= tp,

(5)

f is calculated by using equation (3) for t > t_p, where  tp is the time to surface ponding from rainfall inception.  The relationship between infiltration rate (f) and cummultaive infiltration (F) is:
 

(6)

 or in finite difference form:
 

(7)

or

(8)

where j and j-1 are the time indices,  t_j and t_j-1 are time at the end and beginning of the time step, respectively.

Over a short time period (such as one second), we assume the rainfall is homogeneous and steady in each cell.  Since at time zero the cumulative infiltration is zero, the surface ponding time can be determined at each time step from Equations (3), (5), and (8), by calculating the rainfall rate, surface flow, and cumulative rainfall.  Substituting the infiltration rate f from the previous time step into Equation (8) gives cumulative infiltration at the end of the current time step F_j.  This cumulative infiltration is then substituted into Equation (3) to estimate the new instantaneous infiltration rate.  Ponding occurs when that rate exceeds the sum of the rainfall rate 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 precipitation rate and the water supply from the surface.  Before surface ponding, the infiltration rate is determined by Equation (3) and by the total surface water that can be supplied for infiltration.

In areas where permafrost is near the surface, the cumulative infiltration is limited by the depth of the active soil layer.  The tundra watershed studied by Hinzman and Kane (1992) had continuous permafrost and was found to have an average summer maximum active layer depth of 500 mm with a range of 250-1000 mm.  Consistent with other reports (e.g, Viereck et al., 1986), depth was related to slope, aspect, and drainage characteristics of soil.  Viereck et al. (1986) measured active layer depths in upland black spruce forests near Fairbanks, Alaska, which ranged between 220-510 mm, and between 160-550 mm in lowland black spruce forests.  As the depth of the organic layer increases, soil is insulated and the depth to permafrost is less.  Viereck et al. (1986) found organic layer depths between 170-380 mm and 140-250 mm in upland and lowland black spruce forests, respectively.

We assumed that permafrost varied from the surface to depths of 1000 mm over our site.  The site was initially stratified into floodplain, lowland, and upland regions and permafrost depth was estimated using an empirical relationship following the conceptual model of Viereck et al. (1983). We used aspect as the first determinant of permafrost depth, where north = 160 mm, south  = 1000 mm, and east and west were set at 580 mm.  A linear stretch was used to interpolate permafrost from aspect.  For the upland region, these potential depths were decremented by 10% if the slope was shallow (<5%), and again by 70% if the soil type had poor drainage characteristics and by 10% more if canopy cover was high (>0.27 NDVI).  This allowed the estimated soil active layer in the montane region to range between 40-1000 mm.  The lowland zone used a different model because all soils were classed as poorly drained.  In this zone, the minimum depth in the north aspect was 200 mm and 700 mm in the east and west aspects.  Low slopes decremented the depth by 70% and high canopy cover by an additional 10%.  For the floodplain, all sites had near zero slopes and undefined aspects and all soils were poorly drained.  We used vegetation type and canopy density to estimate depth to permafrost within this zone.  The broadleaf stands (aspen, balsalm poplar (BP), and paper birch(PB)) were given a maximum soil depth of 1000 mm which was decremented by 25% if cover was high.  Black spruce (BS) sites had the least depth with a base of 200 mm at low cover and 160 mm at high cover; white spruce (WS) was assigned a maximum depth of 750 mm under low cover and 500 mm under high cover.   Mixed (M) hardwood and conifer forests were given intermediate values using a linear stretch.  The validity of these assumptions was evaluated by examination of the spatial variation within the site and by randomly extracting pixels and comparing estimated depth and site characteristics with Figure 2 and Table 1 in Viereck et al., 1983.  The predicted soil active layer depths for 33 randomly selected pixels are shown in Table 3.

The soil column becomes saturated when cumulative infiltration is greater than or equal to f-q_i times the soil active layer depth.  Next, the infiltration rate is limited by saturated hydraulic conductivity.  The water naturally steps from the impermeable layer depth to the surface when the wetting front reaches the active layer.  When the subsurface is saturated, exfiltration may occur due to the subsurface flow.  The surface-subsurface flow system is dominated by changes in the hydraulic head and changes in storage.  The initial soil moisture was determined by assuming that the previous summer storm had saturated the soil column and the site had drained for a week prior to the start of the modeled storm.  The initial water content was determined by assuming that the soils were approximately at field capacity.  Because no physical or chemical analyses were available for our soils we examined other soils within Alaska having similar textures and site descriptions.  All five SCS soils were located in the south coast and none was from interior soils (Reiger, 1979).   The mean field capacity for these silty loam soils was 50% +/- 2% volumetric soil moisture, so we assumed this value for field capacity.  Porosity, soil suction head, and saturated hydrologic conductivity were from Rawls and Bakensiek (1983).
 

Subsurface Flow

Boundary and Initial Condition

For determining the boundary condition for each edge and internal cell, we assume that the boundary condition of the interior cell is that determined in the previous time step. Furthermore, storage in a cell is only affected by its four adjacent neighbor cells (in the four cardinal directions), and the flow system is solved explicitly.

Famiglietti and Wood (1994b) suggested initializing the hydrologic model by assuming that the local soil profile state is at gravitational equilibrium which is the assumption we used.  The soil moisture state was initialized by assuming that the soils had drained from saturation for one week prior to the rainfall event.  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 except where a free surface water table existed.  The one layer soil yields step changes in h when the wetting fron advances to the impremeable layer and the new water table depth is the soil surface.  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.
 

Numerical Solution Scheme

Knowing the discharges along the edges of the four adjacent cells from the previous time step, change in groundwater storage in the cell of interest is:
 

(12)

If inflow exceeds outflow, the water table rises or if more water leaves than enters, it falls.  Similarly, change in surface storage is:
 

(13)

Again, if  inflow exceeds outflow, ponding depth increases or it falls if more water leaves than enters.  In both cases, the new depths are used to estimate surface and groundwater flow for the next time step and the process repeats through time.

In the previous discussion, surface and subsurface flow are considered independently because the wetting front is above the water table.  The downward moving wetting front, predicted by the Green and Ampt equation, is initiated at the beginning of infiltration but will change the water table elevation only after the wetting front reaches the water table.  When the wetting front meets the water table, the column is saturated and the water table depth shifts to the ground surface or to the free water surface in the case of ponding.  The short time steps used in the numerical scheme ensures stability.  The equation for estimating groundwater is unchanged except that the flow area is restricted to the average active zone:

(14)

If groundwater inflow exceeds outflow, the change in storage is exfiltration, or if more water leaves than enters, the change in storage is infiltration (at the rate of saturated hydraulic conductivity).  In the latter case, if the net loss exceeds saturated storage, the water table drops below the ground surface.

Having just calculated infiltration or exfiltration, the surface water storage equation remains unchanged except that infiltration is estimated using the subsurface mass balance rather than the Green and Ampt equation. Again, if surface inflow exceeds outflow, ponding depth increases or it falls if more water leaves than enters.

Global mass balance error (Err) is calculated following the same procedure and is summed over the entire flow system.  The absolute global mass balance was calculated as:
 

(15)

Where Q_in is the total recharge (volume water) into the flow system during time period Dt, which is the sum of all the fluxes into the system from precipitation, surface flow, and subsurface flow.  Q_out is the total discharge from the flow system over the time period of Dt, which represents the sum of the fluxes out of the system from the surface and subsurface. The term DS represents the change in system storage over the time period Dt.

CASE STUDY

Information Processing and Model Parameterization

Study Site

To evaluate the empirical equation used to adjust spatial  precipitation patterns we compared the annual precipitation records from 23 years for four central Alaska weather stations (Circle City, Fort Yukon, Healey 2 NW, and McKinley Park) against elevation.   The regression equation was predicted annual precipitation (mm) = 0.3723 (elevation, m) + 167.14 (r² = 0.94).   Since this is was close to the 33 mm/100 m elevation increase reported by Zhang et al. (1987), we adopted their equation.
 

Data

Vegetation Type Distribution, Density, and Land Cover Classification

The second step was to generate a classification land use / vegetation map from the SPOT satellite image data using a supervised maximum likelihood classification to define roads, buildings, water, and other site characteristics that were not defined in the topographically based vegetation classifier.  Much of the landscape is found in early successional stages, particularly since recent fires have removed >80% of the mature forest during the 1980’s (Alaska Fire Service).  Therefore, a map of “actual” vegetation is needed. This was done by adjusting the original topographic classification using a SPOT satellite estimate of green plant cover, termed the Normalized Difference Vegetation Index (NDVI), in order to separate late seral conifer forests from early succession shrub and hardwood dominated forests.  The NDVI is a normalized ratio of red and near-infrared reflectance ((band 3-band 2)/(band 3 + band 2)).  It was used to modify the topographically determined potential vegetation type into the actual type.  Successional sequences followed the conceptual model by Van Cleve and Viereck (1981) that predicts interior Alaskan forest succession through the following stages: shrub and/or herbaceous, followed by deciduous broadleaf forest, with either Black Spruce or White Spruce forests as the climax vegetation type. After the vegetation map was created, the surface roughness was estimated based on the study by Chow (1959) and Barnes (1967).  Density was estimated for each cell by binning the NDVI values into low, medium and high density classes.
 

Soil Distribution Map and Soil Attributes

The soil active layer thickness was generated using vegetation and soil type distributions based on field data from other sites.  Hinzman et al. (1991) found that the typical active layer thickness at Imnavait Creek (4.0° north of our study site) at the end of summer was about 500 mm and ranged from 250 to 1000 mm.  The variation in depth was due primarily to topographic and drainage characteristics, as well as vegetation and soil types. Viereck et al. (1983) and Van Cleve et al. (1986) examined the permafrost depth at Bonanza Creek LTER and defined a relationship between vegetation type and stand age to depth of permafrost.  Bonan (1989, 1992) used similar relationships and climate data to define depth of thaw for forest stands at Fairbanks.  These data give a similar range of active layer depths as Hinzman et al. (1991).  We used the relationships described in Viereck et al. (1983), Van Cleve et al. (1986) and Bonan (1989) to empirically define active layer depths which would coincide with the mid-August storm event (Kane et al., 1991).  Thermal changes in the soil were not considered in the model, so ice lenses inside the active layer are treated as components of the soil.  Table 3 shows the site characteristics and estimated active layer depths for 33 randomly chosen cells from the three elevation zones.
 

Topographic Derived Variables

Simulation and Visualization

The spatially continuous output results can be visually represented in image form and the temporally continuous cell-based output in graphic form.  The full representation of the data output is only realized in a motion picture format.  The entire model output from this study is 289 MB.  Due to the large amount of data, secondary decisions on presentation, validation, and synthesis of the model results were essential.  In total, 10 output layers were computed (from 12 input layers) over 980 time steps for the study area (6,279,000 cells), which required about one  hour to compute on a DEC 3000 workstation.  The time step length was varied with changes in hydrologic processes, from seconds at the early stage of a rainfall event to minute and hour intervals for simulating delayed events like subsurface flow.

The output results from the model are presented in both spatial and temporal forms.  Selected time steps and images were chosen to convey information about the hydrologic processes that characterize this arctic region.  Rationale for the data presentation is described in the next section.

RESULTS AND DISCUSSION

Results from the hydrologic model are presented in Figures 5 to 12 to demonstrate the complex spatial and temporal patterns and their evolution with time.  The study region was divided into flood plain, lowlands (or river terrace), and mountains because of their different precipitation zones and vegetation patterns.  Four cells within the scene were selected for comparison.  P1 is on the flood plain, P2 on the lowland or river terrace, P3 on the north side of the mountain (leeward in this simulation), and P4 on south side of the mountain (windward).  Each point was randomly selected from these topographic categories to represent contrasting hydrologic regimes.   Elevation, terrain features (aspect, slope), soil types, and land cover types of each site (P1, P2, P3, P4) are given in Table 2.

Precipitation increases with elevation and is highest in the south facing direction (Figure 5 and Equation 1).  South facing slopes in the mountains receive significantly greater rainfall than north facing slopes in the lowland areas. For example, solving for the precipitation at P4, which is located at an elevation of 944.9 m and with an aspect of 159^o, we set c=1, and a and b are constants equal to 0.0016 m^-1 and 0.001 degree^-1, respectively.  At one hour the rainfall rate at the weather station is 32 mm/hr so the rainfall rate at P4 is 74 mm/hr.  The rate at P2 on the lowland terraces is higher than P1, corresponding to its higher elevation.  At P3, the rainfall rate is greater than the rainfall rate at P2 and it is highest at P4.  At time = 7 minutes, the rainfall rate is low and ranges from 5 mm/hr in the flood plain site (P1) to 10 mm/hr on the mountain site (P4). As estimated using Equation 1, the rainfall rate increases 2 mm/hr between P1 and P3.  At 30 minutes, the rainfall increases from 16 mm/hr at P1 to 33 mm/hr at P4.  At 60 minutes, the rainfall rate increases in lowland and mountain sites to 32 mm/hr and 74 mm/hr respectively.  After one-half hour, the instantaneous rainfall rates range from 14 mm/hr in lowland  sites to 31 mm/hr on the windward slopes of the mountain sites.

The rainfall patterns vary with time which are shown for the same four locations in Figure 6 (a).  For this typical but large regional summer storm, the maximum rainfall occurred approximately 60 minutes after storm initiation.  The precipitation rate and duration of the August 14, 1987 convective storm is similar to the largest summer storms recorded at the Circle City, Alaska station over the past thirty-five years.  We assumed that all precipitation occurred in two hours.

Infiltration rate at 7, 30, 60, and 90 minutes varies with elevation, soil moisture, aspect, and vegetation cover (Figure 7). Prior to ponding infiltration rate is limited by the rate of water input from precipitation.  At time = 7 minutes, the rainfall rate is low and the soil is relatively dry and the potential infiltration rate predicted by the Green and Ampt equation exceeds the precipitation rate. On the flood plain and on the terraces the infiltration rate is less then 5 mm/hr.  However, the infiltration rate can reach 10 mm/hr in the mountains because precipitation is higher.   In lake beds or river/stream beds, the infiltration rate is less than 2.5 mm/hr because the infiltration rate is limited by the soil moisture content, which is initially affected by water table depth, and later by the lateral subsurface flow.  We followed Famiglietti and Wood (1994) by assuming that the local soil profile is initialized at gravitational equilibrium.  Surfaces covered by lakes or streams may have exfiltration because the underlying soil is saturated.  After 30 minutes, the precipitation rate has increased to nearly maximum, hence, the infiltration rate increases to 5 mm/hr on the lowlands and to 35 mm/hr in the mountains.  After 60 minutes, the infiltration rate is lower in the flood plain (10 mm/hr) than in the mountain region (20 to 35 mm/hr) because most cells in flood plain are saturated.  In the upper mountains, some cells are already saturated or have a deeper wetting front where the initial soil water content and precipitation rate are higher. The spatial infiltration pattern at this time clearly shows that infiltration is related to rainfall rate and changes in soil moisture.  At time = 90 minutes, the soil is saturated in most regions and the infiltration rate reaches an equilibrium status, and is limited only by the subsurface flow velocity and the saturated vertical hydraulic conductivity.  At this time, the spatial pattern of infiltration rates is soil controlled.  The infiltration rate is less than 10 mm/hr for most cells but it can reach 30 mm/hr in regions where soils have a deeper water table and low soil water content.

As mentioned above and shown in theory, Green-Ampt infiltration rate is dominated by soil physical properties. The coarser the soil, the higher the infiltration rate. In this study region, the soil type varies from silt loam in flood plains to sand loam in the mountains.  However, due to lack of site-specific physical data, we assumed all six soil types had the same physical characteristics.  So in this example, the infiltration rate is controlled by rainfall during the early stages of a storm event.  The first site, P1, has the lowest infiltration rate among the four sites due its lower rainfall rate.  Site P2 has a slightly higher rainfall rate than P1, and P4 which has the highest elevation and is south facing, has the highest rainfall rate and highest infiltration rate.   Because elevation at P3 is between P1 and P4, its infiltration rate is intermediate.  The significant influence of rainfall rate on infiltration rate and subsequent changes in infiltration with time are shown in Figure 6 b.  In the early stages of the rainfall event, the predicted infiltration rate from the Green-Ampt equation is greater than the rainfall rate.  Since the actual infiltration rate is controlled by rainfall rate the infiltration rate increases with time.  Later the rainfall rate exceeds the infiltration rate due to soil saturation and the predicted infiltration rate is less than the rainfall rate.  Infiltration generally decreases with time.  After the rainfall event stops, the infiltration processes are dominated by subsurface flow and surface water ponding.

The spatial distribution of cumulative infiltration is similar to that of infiltration rate (Figure 7) and not shown here.  The infiltration rate and cumulative infiltration are important factors in understanding pollution and irrigation 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 8 shows spatial patterns in surface ponding at 7 minutes and at 2 hours, 5 hours, and 12 hours.  Times beyond the two-hour rainfall event were selected to illustrate the persistence in surface flow.  From these images, one observes the strong dependence of surface ponding on terrain features.  At 7 minutes, surface ponding depth is less than 2.5 mm in all cells.  At 2 hours, the rainfall event stopped.  The surface water ponding depth in lowlands is between 20 and 40 mm but for the mountains the ponding depth can reach 400 mm in the stream channels. Runoff from the hill slope region fills low lying areas.  At 5 hours, the surface ponding depth in sites having steep slopes has decreased, but the depth in stream channels continues to increase.  Later, at 12 hours, most surface ponding is restricted to stream channels. Changes in surface water ponding depth over time are shown in (Figure 6 c).  Surface ponding at  P4 is shallow and only occurs during the earliest times in the simulation.  In contrast, the pond depth continues to increase in lower elevation areas as more water is transported by overland flow from the mountains. P3 has the highest ponding depth and reaches a maximum after 4 hours while the change in ponding depth with time at P4 is more related to the time required for flow from its contributing area.  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 9 shows the spatial patterns in surface flow velocitys at 7 minutes, 2 hours, 5 hours, and 12 hours.  At 7 minutes, the surface flow velocitys are low (0.005  m/s) due to the limited total rainfall.  At 2 hours, higher overland flow occurs in most cells.  In the mountain regions, because of the steep terrain and high rainfall rates the surface flow velocitys are higher than in the lowlands. The surface flow velocity is equal to or greater than 0.1 m/s on all sites except the mountain ridges which have partially drained. At P1 in the flood plain, the flow velocity is low (0.01 m/s) due to its low elevation and slope and the limited surface water ponding depth.  At P2, the surface flow velocity reaches 0.01 m/s, because the ponding depth is low even though the slope is steeper than at P1. The same reasons explain the higher surface flow velocitys at P3 and P4. However, at 5 hours overland flows have largely stopped in mountainous areas except in the stream channels while the flow velocitys are still increasing in the lowlands.  At 12 hours, surface flows have been mostly confined to the stream channels and flood plains.  The stream distribution and contributing area are clearly shown in the spatial flow patterns.

From Figure 6 d, we see that the peak flow and the flow velocitys vary over time at the four sites.  Early in the rainfall event. P1 had a peak surface flow velocity of approximately 0.02 m/s which is the lowest peak flow and occurs at 21 hours.  The peak flow at P2 is 0.03 m/s and is reached quickly;  high flow velocitys last the briefest time among the four sites.  P3, a site with high vegetation cover, reaches peak surface flow early but the flow is distributed more evenly over a longer time interval.  Lastly, P4 was located close to a stream channel in the mountains and as expected, the peak surface flow occurred at the time the rainfall stopped and the rate is higher than at P3.

The water table depth shown in Figure 10 demonstrates the variations in subsurface water resources with time. We use the ground surface as the reference data for water table depth.  If the soil is unsaturated to the permafrost, water table depth is equal to the depth of the active layer.  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 minutes, little precipitation has fallen and the water table is deep (common in the mountains) except where cells are already occupied by surface water bodies or are flat regions where the soil is saturated. At 2 hours, all cells become saturated and the water table approaches the ground surface.  By 12 hours, most cells are saturated.  On the upper mountain regions, due to higher rates of subsurface flow, the water table falls below the ground surface.  This pattern is more prominent after 54 hours. Figure 6 e shows how the water table changes with time. As time increases, the water table increases by the water contributed from infiltration and subsurface flow  and drops by the water  contributed to subsurface flow. P1 has a shallow soil active layer depth and becomes saturated after 7 minutes.  P2 has the same pattern as P1, however, the impermeable layer at P3 occurs at much greater depth.  Water table changes at P4 are similar to those at P3 except for the greater active layer depth.

Subsurface flow velocity is important for contamination studies as well as characterizing water resources distribution.  Figure 11 shows the pattern of subsurface flow at 7 minutes, 2, 12, and 54 hours.  At 7 minutes, most cells in the lowland and flood plains have subsurface flow velocitys less than 0.05 m/day. Where saturated soils exist and slopes are steep, subsurface flows reaches 1.8 m/day. Several artifacts are observed in Figure 11 as yellow and red class geometric shapes in the flood plain and lowland regions.  These facies are produced by the estimate of active layer depth and are attributable to the resolution of the DEM in sites having low slopes.  At 2 hours, the subsurface flows have increased over most of the area and higher subsurface flow can be found in the mountain region on steep slopes. The mountain region receives more rainfall than the lowlands and the flood plain, and therefore, the water table is rises more quickly.  At 12 hours, subsurface flow occurs in most cells, even in the upper regions of the mountain. Subsurface flow is low only in the lowest regions where the hydraulic gradient is low. The same pattern occurs at 54 hours.   Figure 6 f  shows that the subsurface flow velocity varies with time at the locations P1, P2, P3, and P4.  The very different subsurface flow patterns observed among these sites were largely affected by topography, especially in sites with steep slopes.
 

Mass Balance Error Estimate

SUMMARY

The specific results in this study have not been rigorously tested because of a lack of sufficiently detailed spatially explicit ground observational data in the study region.  Such detailed data would be difficult to obtain nearly anywhere, except in gauged catchments, and are particularly unavailable in relatively pristine areas like this watershed where surface access is nearly impossible due to bogs and thermokarst terrain.  Nonetheless, because better understanding of hydrologic processes in arctic regions is of critical importance for many ecological and climate issues, spatially explicit models that utilize available environmental data bases are needed. Modeling of large scale processes must be validated through a variety of indirect techniques including simulation and testing of the individual submodel components, evaluation of the results as a system (i.e., for the system as a whole to be correct, all components must be within realistic ranges).  Evaluation of the spatial image maps of each of the parameter outputs and through the temporal hydrographs from the four sites indicated that all sites produced spatial and temporal patterns consistent with theoretical expectations.

Because the model currently lacks an evapotranspiration estimation function, it can only be applied to storm event simulations. We plan to add evapotranspiration and surface energy budget functions to the model (e.g., FOREST-BGC (Running and Coughlan, 1988) or SiB (Sellers et al., 1986)), so that it can be used to study multi-event or seasonal simulations.  Inter-seasonal simulations also require incorporation of snow pack and snow melt processes.  If the model were applied to other sites to estimate ET, a multiple-layer soil characterization (e.g., that of Wigomosta et al., 1994) would improve the results.  This also assumes that the vegetation map can be used to estimate rooting depth, which is needed for calculating soil moisture redistribution, evaporation, and transpiration.  The accuracy of model initialization, particularly soil moisture, could be improved by running the model for an extended period of time (one to two years) with true climate data as input.
 

Acknowledgments

REFERENCES

ARC /Info User's Guide.  Environmental Systems Research Institute, Inc. Version 6.1. Apr., 1992.

Band, L. E.,  Effect of land surface representation on forest water and carbon budgets, J. Hydrol., 150, 749-772, 1993.

Band, L. E., Patterson, P., Nemani, R., and Running, S. W., Forest ecosystem processes at the watershed scale: Incorporating hillslope hydrology, Agri. For. Meteorol., 63, 93-126, 1993.

Barnes, H. H., Roughness Characteristics of Natural Channels,  U.S. Geological Survey Water Supply Paper 1899, 1967.

Batelaan, O., De Smedt, F., Otero Valle, M. N.,  and Huybrechts, W.,  Development and application of a groundwater model integrated in the GIS GRASS,  HydroGIS 93: Application of Geographic Information System in Hydrology and Water Resources. Kovar, K., and Nachtnebel, H.P. (eds.), Int. Assn. Hydro. Sci. Publ. (IAHS), Wallingford, Oxfordshire, OX10 8BB, UK, No. 211. pp. 581-589, 1993.

Bergström, S., Development and application of a conceptual runoff model for Scandinavian catchments, Rep. RH07, 118pp., Swedish Meteorol. and Hydrol. Inst., Norrköping, Sweden, 1976.

Beven, K.J. and Kirkby, M.J., A physically based, variable contributing area model of basin hydrology, Hydrol. Sci. Bull., 24, 43-69, 1979.

Bonan, G. B., A computer model of the solar radiation, soil moisture, and soil thermal regimes in boreal forests.  Ecol. Model. 45, 275-306, 1989.

Bonan, G.B., Soil temperature as an ecological factor in boreal forests, A systems analysis of the global boreal forest, Shugart, H.H., Leemans, R., and Bonan, G.B. (eds.), Cambridge University Press, New York, pp. 126-143, 1992

Chow, V. T.,  Open Channel Hydraulics,  McGraw-Hill,  New York,  705 pp., 1959.

Dingman, S.L., Characteristics of summer runoff from a small watershed in central Alaska, Water Resour. Res., 2, 751-754, 1966.

DeVries, J. J., and Hromadka T.V.,  Computer Models for Surface Water Handbook of  Hydrology,  pp. 21.1-21.39,  1993.

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

Famiglietti, J. S., Wood, E.F. , Sivapalan, M., and Thongs, D.J.,  A catchment scale water balance model for FIFE, J. Geophys. Res., 97 (D17) 18, 997-19,007, 1992.

Famiglietti, J. S., and Wood, E.F., Multiscale modeling of spatially variable water and energy balance processes, Water Resour. Res., 30, 3061-3078, 1994a.

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

Famiglietti, J. S., and Wood, E.F., Effects of spatial variability and scale on areally averaged evapotranspiration, Water Resour. Res., 31, 699-712, 1995.

Gao, X., Sorooshian, S.  and Goodrich, D. C.,  Linkage of a GIS to a distributed rainfall-runoff model,  Environmental Modeling with GIS.  Goodchild, M.F., Parks, B.O., Steyaert, L.T., (eds.), Oxford Press, New York, pp. 182-187, 1993.

Goodrich, D. C., Geometric simplification of a distributed rainfall-runoff model over a range of basin scales,  PhD Diss., Dept. of Hydrology and Water Resources, University of Arizona, Tucson, AZ,  361pp, 1990.

Green, W.H., and Ampt G.A., Studies on Soil Physics: 1. Flow of air and water through soil,  J. Agric. Sci.,  4,  1-24,  1911.

Hinzman L.D., and Kane, D.L., Snow hydrology of a headwater arctic basin,  2. Conceptual analysis and computer modeling, Water Resour. Res., 27, 1111-1121, 1991.

Hinzman L.D., and Kane, D.L.,  Potential response of an arctic watershed during a period of global warming, J.  Geophys. Res., 97 (D3), 2811-2820, 1992.

Hinzman L. D., Kane, D.L, Gieck, R.E., and Everett, K.R., Hydrologic and thermal properties of the active layer in the Alaska arctic,  Cold Regions Sci. Technol., 19, 95-110, 1991.

IDL User's Guide.  Research System, Inc.,  V3.5.  June, 1993, Boulder, CO 80303, 1993.

Kaden, S. O.,  GIS in water-related environmental planning and management:  Problems and solutions,  HydroGIS 93: Application of Geographic Information System in Hydrology and Water Resources, Kovar, K., and Nachtnebel, H.P. (eds.), Int. Assn. Hydro. Sci. Publ. (IAHS), Wallingford, Oxfordshire, OX10 8BB, UK, No. 211. pp. 385-397, 1993.

Kane D. L., Hinzman, L. D., Benson, C. S.,  and Everett, K. R.,  Hydrology of Imnavait Creek, an arctic watershed, Holarctic Hydrology, 12, 262-269, 1989.

Kane D. L., Gieck, R.E., and Hinzman, L. D., Evapotranspiration from a small Alaska arctic watershed, Nordic Hydrol., 21, 253-272, 1990.

Kane D.L., Hinzman, D.L., Benson, C.S., and Liston, G.E., Snow hydrology of a headwater arctic basin 1. Physical measurements and process studies, Water Resour. Res., 27, 1099-1109, 1991a.

Kane D.L., Hinzman, D.L., and Zarling, J.P., Thermal response of the active layer to climatic warming in a permafrost environment, Cold Regions Sci. Technol., 19, 111-122, 1991b.

Leaversley, G. H., Lichty, R.W., Troutman, B. M., and Saindon, L.G., Precipitation-runoff modeling system--users manual.  Water Resources Investigations Rep. 83-4238, Washington, D.C., US Geol. Surv., 207pp., 1983.

Maidment, D. R., GIS and hydrological modeling, Environmental Modeling with GIS.  Goodchild, M. F., Parks, B. O., Steyaert, L. T., (eds.), Oxford Press, New York, pp. 147-167, 1993.

Mein, R.G., and Larson C. L.,  Modeling infiltration during a steady rain, Water Resour. Res., 9, 384-394,  1973.

Moore, I. D., Turner, A. K., Wilson, J. P., Jenson, S. K., and Band, L. W., GIS and land-surface-subsurface process modeling, Environmental Modeling with GIS.  Goodchild, M. F., Parks, B. O., Steyaert, L. T., (eds.), Oxford Press, New York, pp. 196-230, 1993.

Nachtnebel, Furst, H. P. J., and Holzmann, H.,  Application of geographical information system to support ground water modeling,  HydroGIS: Application of Geographic Information System in Hydrology and Water Resources Management. Kovar, K., and Nachtnebel, H.P. (eds.), Int. Assn. Hydro. Sci. Publ. (IAHS), Wallingford, Oxfordshire, OX10 8BB, UK, No. 211. pp. 653-664,  1993.

Nemani R., Running, S. W., Band, L. E., and Peterson, D. L., Regional hydroecological simulation system:  An illustration of the integration of ecosystem models in a GIS,  Environmental Modeling with GIS, Goodchild, M.F., Parks, B.O., Steyaert, L.T., (eds.), Oxford Press, New York, pp 296-304,  1993.

Ostendorf  B., and Reynolds, J. F., Relationships between a terrain-based hydrologic model and patch-scale vegetation patterns in an arctic tundra landscape, Landscape Ecol., 8, 229-237, 1993.

PACE Reference Manual.   PCI, Inc.  Version 5.0a.  December, 1991, Richmond Hill, Ontario, Canada, L4B1M5, 1991.

Patric, J.H., and Black, P.E., Potential evapotranspiration and climate in Alaska by Thornthwaite’s classification.  USDA Forest Service Research Paper PNW-71.  Pacific Northwest Forest and Range Experiment Station, Portland, OR, 28p., 1968.

Rawls, W.J., Brakensiek, D.L., and Saxton, K.E., Estimation of soil water properties, Trans. Am. Soc. Agric. Eng., 25, 1316-1320, 1328, 1982.

Rawls, W.J. and Brakensiek, D. L.,  A procedure to predict Green and Ampt infiltration parameters,  Advances in Infiltration.  Am. Soc. Agric. Eng.,  pp. 102-112,  1983.

Rieger, S., Schoephorster, D.B.,  and. Furbush, C. E.,  Exploratory Soil Survey of Alaska,  National Cooperative Soil Survey, USDA,  Soil Conservation Service, February, 1979.

Ross, B. B., Contractor, D. N., and Shanholtz, V.O., A finite element model of overland and channel flow for assessing the hydrologic impact of landuse change, J. Hydrol., 41, 1-30, 1979.

Running, S.W., Nemani, R.R., and Hungerford, R.D.,  Extrapolation of synoptic meteorological data in mountainous terrain, and its use for simulating forest evapotranspiration and photosynthesis.  Can. J. Fort Res., 17, 472-483, 1987.

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

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

SIPS User's Guide, V1.1.  Center for the Study of Earth From Space, Cooperative Institute for Research in Environmental Science.  University of Colorado, Boulder.   March, 1992.

Starosolszky, O. (ed.),  Applied Surface Hydrology.  Water Resources Publications, Littleton, CO 80161,  pp. 40-120, 1987.

Takhtajan, A., Floristic Regions of the World,  University of California Press, Berkeley, CA, 1986.

Thornthwaite, C. W., An approach toward a rational classification of climate, Geograph. Rev., 38, 55-95, 1948.

Ustin, S. L.,  Szeto, L.-H., Xiao, A.-F., Hart, Q.J. and Kasischke, E.S., Vegetation mapping of forested ecosystems in central Alaska, in Int. Geosci. and Remote Sens. Symp., IGARSS-94, Pasadena, CA, August 8-12, 1994.

Van Cleve, K., and Viereck, L.A.,  Forest succession in relation to nutrient cycling in the boreal forest of Alaska.  In:  West, D.C., Shugart, H.H. and Botkin, D.B. (eds), Forest Succession: Concepts and Applications, pp. 185-211, Springer-Verlag,  NY, 1981.

Van Cleve, K., Dyrness, C. T., Viereck, L. A., Fox,  J., Chapin, F. S. III, and Oechel, W.,  Taiga ecosystems in interior Alaska, Biosci., 33, 39-44, 1983.

Van Cleve, K., and Yarie, J., Interaction of temperature, moisture, and soil chemistry in controlling nutirent cycling and ecosystem development in the tiaga of Alaska, Forest Ecosystems in the Alaskan Taiga, A Synthesis of Structure and Function, Van Cleve, K, Chapin, F.S., III., Flanagan, P.W., Viereck, L.A., and Dyrness, C.T. (Eds.), Springer-Verlag, New York, pp. 160-189, 1986.

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

Viereck, L. A., Dyrness, C.T., Van Cleve, K., and Foote, M. J., Vegetation, soils, and forest productivity in selected forest types in interior Alaska, Can J. For. Res. 13, 703-720, 1983.

Viereck, L.A., Van Cleve, K., and Dyrness, C.T., Forest ecosystem distribution in the taiga environment,  Forest Ecosystems in the Alaskan Taiga, A Synthesis of Structure and Function, Van Cleve, K, Chapin, F.S., III., Flanagan, P.W., Viereck, L.A., and Dyrness, C.T. (Eds.), Springer-Verlag, New York, pp. 22-43, 1986.

Vieux, B.E.,  Geographic Information System and non-point sources water quality and quantity modeling,  Hydrol. Process.,  5, 101-113,  1991.

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

Xiao, Q.F., A spatially and temporally continuous surface-subsurface hydrologic model, University of California, Davis,  MS Thesis, 123pp., 1994.

Zhang, J., and Deng, Z.,  The relationship between topography and precipitation in mountain area--Precipitation Study of XinJiang, China,  pp. 160-168. Meteorology Press, Beijing, China, 1987.