David R. Maidment

January 22-26, 1996

Santa Fe, New Mexico

- Developing a Spatial Hydrology Model
- Time Averaged Hydrologic Modeling
- Time Varying Water Balance Models
- Conclusions
- References

The First International Conference on GIS and Environmental Modeling was held in Boulder, Colorado, in September 1991. At that Conference I presented a survey of the state of GIS and hydrologic modeling as the subject then existed (Maidment, 1993). The intent of this paper is to measure the progress of the subject in the four years that have elapsed since the first conference, to present a framework within which GIS hydrology can be viewed as an integrated subject, and to examine how spatial hydrologic models can be created for time averaged and time varying systems. The focus of this paper is on a personal approach to this field rather than a survey of all the relevant applications.

Important progress has been made during the last four years in the availability of comprehensive spatial data sets which support hydrologic modeling. The distribution of low cost or free data sets via Internet and CD-ROM for digital elevation data, soils, land use, and climate data has stimulated the development of procedures for processing the data into useful forms for hydrology, such as those shown on our GIS hydrology home page (http://www.ce.utexas.edu/prof/maidment/gishydro/). It has also stimulated the construction of a number of integrated systems where existing hydrologic models have been connected to spatial data bases resident in GIS, with the SWAT (Soil Water Assessment Tool) of the USDA and the MMS (Modular Modeling System) of the USGS being probably the most comprehensive systems of this type.

A further and more complicated question is to ask how hydrologic modeling can be rethought in the spatial context that GIS provides. In other words, instead of attaching existing models to GIS databases, can new hydrologic models be created that take advantage of the spatial data organizing capabilities of GIS? This question implies a reversal of traditional priorities in hydrologic modeling where the emphasis has always been on the way that physical processes are represented, and the manner in which the parameters are to be obtained for a particular environment plays a relatively minor role. In a spatial hydrology model, the emphasis is first on the digital description of the environment, and then on the formulation of process models which can fit the available data and environmental description.

Hydrology is concerned with study of the motion of the earth's waters through the hydrologic cycle, and the transport of constituents such as sediment and pollutants in the water as it flows. GIS is focused on representing the landscape by means of locationally referenced data describing the character and shape of geographic features. A spatial hydrology model is one which simulates the water flow and transport on a specified region of the earth using GIS data structures. Suppose the boundary of this region is represented by a polygon, such as a river basin boundary or an aquifer boundary. Because both vertical and horizontal water flow can be taking place within the region, hydrologic processes need to be defined over a volume of space rather than an area. Such a volume can be constructed by projecting vertically the lines making up the polygon boundary into the atmosphere and into the earth, and closing the top and bottom of the volume by horizontal planes. Such a volume is called a *control volume* in fluid mechanics and the surface which surrounds it is called the *control surface*.

Building a hydrologic model involves writing equations that relate the rates of change of water properties within the control volume to flow of those properties across the control surface. For example, a simple soil water balance model for a control volume drawn around a block of soil is:

S(t+1) = S(t)+P(t)-E(t)-Q(t)

in which S(t) represents the amount of soil moisture stored at the beginning of the time interval t, S(t+1) the storage at the end of that interval, and the flow across the control surface during the interval consists of precipitation P(t), evaporation, E(t), and soil moisture surplus, Q(t), which supplies streamflow and groundwater recharge. Solving this equation requires dealing with time series of the four variables: S, P, E, Q, and possibly of other variables related to them.

It is implicit in constructing a spatial hydrology model that the properties of the system will be spatially variable, so a time series for each of the variables just described must be generated for each soil unit in the domain of analysis. Suppose that there are L spatial units, and that analysis on a single unit requires the definition of M variables for each of N time periods. The number of values to be determined is given by the product LMN, a number that can easily explode beyond the capabilities of available computer memory and reasonable computation times. Indeed, it appears that if the product LMN exceeds a limit of approximately 1 million, solution of the model will be computationally difficult. It follows that in constructing a GIS hydrology model the first task is to determine what variables will be calculated on how many spatial units for a defined number of time periods. Construction of complex models must proceed by partitioning the total problem into a series of submodels that interact with a common database. It is the capacity of GIS for rigorously defining this database which makes possible complex models connecting various parts of the hydrologic cycle within a particular region. One begins with a GIS description of the region and then by modeling adds additional detail to the regional description concerning water flow and constituent transport.

In considering these questions, a fundamental distinction arises in the treatment of spatial units depending upon whether raster or vector data structures are used. In a raster data structure, a fine mesh of cells is laid over the landscape and all calculations are done for each cell. This approach, based on using digital elevation model cells as the spatial units, is very useful for certain kinds of hydrologic analysis. But the number of DEM cells within typical analysis regions is usually very large, typically 10,000 to 1 million, so the number of time periods that can analyzed is usually relatively small. Indeed, one may eliminate time as a dimension of the problem by working with mean annual values of all the variables, and the subject of raster based mean annual flow and transport models is presented at further length later in this paper.

If time dynamics are to be considered, it is usually necessary to employ a vector data structure based on related points, lines and polygons. To describe the values on L spatial units of M computed variables in N time periods requires in concept a 3-D data structure but this can be reduced to a set of 2-D data structures in the following manner. The feature attribute table of the GIS coverage defines the geographic properties of the L spatial units and gives each a unique identifying number. This table has geographic attributes as its columns or fields for which the values in each spatial unit are displayed in rows. The values of a particular computed variable, such as soil moisture storage, S, can be defined by means of a related *time table*. The conventional method of constructing such a table is to define a new field for each time interval and keep the rows for the sequence of spatial units, but the advent of object-oriented GIS programming languages, such as the Avenue language in Arcview, has made possible the reverse arrangement, namely the use of time as the index on the rows of the table for which there is a new field for each spatial unit. The item name in this field is the feature identification number of the spatial unit attached to an arbitrary prefix such as SU to mean spatial unit. Thus, in a time table for soil moisture storage, the value of soil moisture storage in time period 154 on spatial unit 30 may be found in row 154 of the table in field SU30.

This arrangement of time vertically and space horizontally makes it easy to see by reading down a column the temporal sequence of values of a variable in a particular spatial unit. The key point in this description is that using object oriented data structures, a value in one table (the feature identification number) can be related to a field in another table instead of to a single value or a set of such values. For hydrologic modeling, the smooth treatment of time variation within GIS is a critical problem that has strongly limited what could be done in the past. It appears that this limitation has now been lifted and the models of spatially distributed and time varying systems can be readily constructed within GIS rather than simply using the GIS as a repository for spatial data feeding an external time varying model.

In Table 1 is presented a 10-step plan by which preparing a GIS hydrology model can be broken down into component parts. The first five of these steps deal with defining the framework of the model in space, and time, and in preparing the environmental description, which may include representation of the land surface terrain, soils and land cover, subsurface hydrogeology, and hydrologic data such as precipitation, streamflow and constituent concentrations. The second five steps deal with simulating the water balance of spatial units, the flow of water and transport of constituents between units, the effect of water utilization structures such as dams and pumping systems, and finally, with the presentation of the study results.

**Study design:**Objectives and scope of study; spatial and time domain; process models needed, variables to be computed.**Terrain analysis:**Deriving a watershed and stream network layout from digital elevation data and mapped streams.**Land surface:**Describing soils, land cover, land use, cities, and roads.**Subsurface:**Hydrogeologic description of aquifers**Hydrologic data:**Locating point gages, attaching time series and their average values, interpolating point climatic data onto grids.**Soil water balance:**Partitioning precipitation into evaporation, groundwater recharge and surface runoff; partitioning of chemicals applied to the land surface.**Water flow**Movement of water through the landscape in streams and aquifers. Computing streamflow and groundwater flow rates.**Constituent transport:**Transport of sediment and contaminants in water as it flows. Computing concentrations and loadings.**Impact of water utilization:**Locating reservoirs, water withdrawals and discharges in rivers, and aquifer pumping. Their effects on water flow and constitutent transport.**Presentation of results:**Developing visual and tabular presentation of the study results. Use of Internet and CD-ROM to transmit results.

Table 1. A ten-step procedure for a GIS hydrology study.

In our GIS hydrology home page is presented sources of spatial data on Internet that can be used in preparing the environmental description. The advent of STATSGO soil data base and the Internet access to the RF1 river reach file, the USGS Hydrologic Unit Code for watershed boundaries and the USGS Land Use and Land Cover data are all a great help in constructing hydrologic models.

Probably the principal advance in GIS hydrology modeling that has occurred during the past several years has been the widespread availability of digital elevation data via Internet and CD-ROM and advances in the methods of processing them. DEM data are available in the United States at 30m and 3" cell size, suitable for delineating watersheds and stream networks within urban areas or in small river basins, at 15" cell size which is suitable for regional studies of the size of Texas, to 30" cell size which is appropriate for continental scale studies to 5' cell size covering the earth. The 30" DEM's for Africa and other regions being constructed by Sue Jenson and her co-workers at the USGS using the Digital Chart of the World are a significant contribution.

Delineation of watersheds from DEM data has become standardized on the eight-direction pour point model in which each cell is connected to one of its eight neighbor cells (four on the principal axes, four on the diagonals) according to the direction of steepest descent. Given an elevation grid, a grid of flow directions is constructed, and from this is derived a grid of flow accumulation, counting the number of cells upstream of a given cell. Streams are identified as lines of cells whose flow accumulation exceeds a specified number of cells and thus a specified upstream drainage area. Watersheds are identified as the set of all cells draining through a given cell. The "thousand-million" rule is a rough guide in this activity - take the area of the region to be analyzed and divide it by one million to give the appropriate cell size to use; multiply the cell size chosen by one thousand and that is the minimum drainage area of watersheds that should be delineated from this DEM.

Table 2 shows the cell sizes of digital elevation data and their typical range of application. In this table, a typical watershed area contains 5000 DEM cells and the region of application contains 200 of these watersheds.

Geographic Linear Cell Watershed Area Region Area Typical Cell Size Size (km2) (km2) Application

1" 30m 5 1000 Urban watersheds 3" 90m 40 8000 Rural watersheds 15" 460m 1000 200,000 River basins 30" 930m 4000 900,000 Nations 3' 5.6 km 150,000 30,000,000 Continents 5' 9.3 km 400,000 90,000,000 Global

Table 2. Digital elevation cell sizes and their scope of application

A standardized way of delineating watersheds and stream networks is the following: construct the stream network in the standard way by choosing a drainage area threshold; divide the stream network so created into individual stream links; find the outlet cell at the lower end of each link; delineate the watershed for each of these outlet cells. By changing the threshold drainage area, subwatersheds can be delineated within watersheds in a nicely scaled manner. This algorithm has a critical property - there is one and only one stream for each watershed, which makes modeling of the flow of water from watersheds to the streams within them feasible within GIS. Also, there is a one to one relation between the raster representation of the landscape and the vector features of streams and watersheds derived from it.

Several variations on this standard algorithm are useful. First, the DEM delineated streams are usually close to but do not coincide with mapped streams. Sometimes, critical errors occur when a portion of the upstream part of a basin drains into the wrong downstream river. To overcome these errors, the mapped streams can be converted into a grid and "burned in" to the DEM by artificially raising the elevation of the off-stream cells. This technique requires editing of the stream network to eliminate stray streams and loops that would confuse the delineation process and it results in some distortion of the watershed boundaries in areas where the DEM and mapped streams are not completely consistent, but it has the great advantage that the DEM delineated streams match the mapped streams exactly. After all, it is the stream network which is really the critical item in landscape delineation because it is the stream that carries the water. This "burn in" technique is especially useful in coastal zones with very flat terrain and other locations where drainage is directed through constructed channels. It also helps to ensure that gaging stations and other features are precisely located on the stream. Other variations on the standard delineation technique include the identification of zones of interior drainage and the subdivision of long streams by placing outlet cells along the streams at arbitrarily defined intervals. Arc Macro Language Scripts for watershed delineation are available on our GIS hydrology home page

The determination of flow at ungauged locations is a common problem in hydrology. A simple approach to this problem is to eliminate time as a dimension by restricting the computation to mean annual flows. The analysis can then be constructed by using the cells of a DEM grid as the computational units. One begins with a mean annual precipitation grid over the landscape, which for the United States has been constructed by Daly et al. (1994) and for Africa by Hutchinson et al. (1995), both using approximately 3' cells. The precipitation for each DEM cell is determined from the climate grid. The watersheds of each of the stream gauging stations in the region are delineated and the mean annual precipitation, P, for the drainage areas determined. The longest streamflow record in the basin is used as an anchor record, a long period of analysis is chosen (such as 1961-1990), and the mean annual flow per unit of drainage area, Q, is determined for each gage. If some of the gages have incomplete records, the long term estimate of the mean annual flow can be found by: long term flow at a sample gage = long term flow at the anchor gage x (flow at sample gage / flow at anchor gage) where the ratio in parentheses is constructed using the means of the common period of record at the two gages. A graph is plotted of Q versus P, and an equation fitted to the relation:

Q = c P

where c is a runoff coefficient, which is a nonlinear function of the precipitation. In dry areas, the greater is the precipitation, the greater is the percentage of the precipitation which becomes runoff. By multiplying the mean annual precipitation grid by this runoff coefficient, a mean annual runoff per unit area can be determined for each DEM cell. This quantity can be used as a weight and a weighted flow accumulation performed in the same manner as the regular flow accumulation is done when constructing the watershed boundaries. The weighted flow accumulation of each DEM cell, when multiplied by the cell area, gives the mean annual flow for each cell. Thus a mean annual flow map can be derived with estimates of the flow at every stream location in the landscape. This is a very simplified method of hydrologic analysis but one that is faithful to the gauged data in the region and can be applied to large regions in a consistent manner.

An extension of the mean annual flow estimate just derived can be used to estimate non-point sources of pollution loading to streams. Such sources include pollutants from agricultural areas and from urban runoff from areas such as roads and parking lots. Point pollutant sources are those associated with a particular outlet location, such as the outlet of a wastewater treatment plant. A standard assumption in treating non-point source pollution is to relate the expected concentration of pollutants in runoff to the land use in the drainage basin. By taking a land use map of the area, and using a look-up table to connect land use to expected pollutant concentration, a map of expected concentration, C, can be determined and the expected pollutant concentration in runoff from each DEM cell calculated. By taking the product of the concentration and the flow rate per unit area, a pollutant loading rate per unit area, L, can be calculated for each cell using the relationship:

L = CQ

Computing the weighted flow accumulation of this loading onto downstream cells and multiplying by the cell area gives the mean annual pollutant loading in each cell. Nice maps of expected pollutant loadings can thus be created. By dividing this mean annual loading by the mean annual flow, a grid of expected pollutant concentration in each cell is derived, which is derived from the weighted average of all the contributions of flow and loading from upstream. Thus is created an expected concentration map, which of course is really meaningful only in the stream cells.

Pollutant concentration is normally sampled at a number of locations in the basin. By making the assumption that the concentration sampled each time at a particular place is drawn from the same water (i.e. that the data are statistically stationary in time), an average observed concentration can be computed for each sample point. The observed average concentrations can be compared with the expected mean annual concentrations computed from the flow accumulation model by drawing a circle at each sample point having an area proportional to the number of samples, and using a consistent color coding scheme for the expected and observed average concentrations. This analysis quickly shows where pollution is at or near expected levels, where point sources are raising the concentration well above the levels expected from non-point sources alone and provides a spatial picture of the variation of observed pollutant levels. Like the mean annual flow analysis, this grid-based non-point source pollution assessment method is simplified but it makes use of the observed hydrologic data and the GIS data normally available in a region in a reasonable way and yields useful results.

A *water balance model* is a representation of the mass balance of water within a particular control volume. It is a physical statement of the law of conservation of mass which states that matter cannot be created or destroyed. As a result, the rate of change of storage of water within the control volume is equal to the difference between its rates of inflow and outflow across the control surface. One may distinguish in constructing a spatial hydrology model between the surface defining the outer boundary of the study region, and the surface defining the boundary of the spatial units within that region. A spatially distributed water balance model applies the law of conservation of mass to describe the mass balance within each spatial unit, and to this must be coupled a momentum equation (such as Darcy's law for groundwater flow) which defines how quickly water can move between units. Different sizes and shapes of spatial units are needed to deal with the different phases of the hydrologic cycle.

Most water that falls on the land surface is derived from oceanic evaporation carried inland by atmospheric circulation, so it is appropriate to begin the study of hydrology by examining the motion of atmospheric water. The most useful way of doing this in a GIS context is to use the results of GCM modeling, where the acronym GCM means here General Circulation Model (this was the original meaning of this acronym before the more popular Global Climate Model came into vogue). In the United States, the National Meteorological Center in Maryland maintains a global GCM in continuous operation for numerical weather forecasting, which is updated each 12 hours with data from atmospheric soundings obtained from a global network of balloon-borne sensors released from weather stations, called the Global Data Assimilation System. The condition of the atmosphere (temperature, density, wind velocity, air pressure, moisture content) is calculated on a geographic grid of 2 degree cells covering the earth, using a very short time interval of the order of a few minutes, for a time horizon of a few days ahead. Each 12 hours, the forecasts are updated with new observations, and the process repeated. Summary statistics are available from the National Center for Atmospheric Research in Boulder, Colorado. A similar service is provided in Europe by the European Center for Medium Range Weather Forecasting in Reading, England.

An example of the atmospheric moisture flow over North America is provided in our GIS hydrology home page, with atmospheric moisture flow being given by the product of the wind velocity and atmospheric moisture content, vertically integrated over the height of the atmosphere, spatially averaged over the 2 degree grid, and time averaged over a month, for the months June 1991 to August 1993. An animation of these data provided at this home page site, shows that the bulk of the water precipitating over the central and eastern regions of the United States is derived from evaporation in the eastern Caribbean Sea, carried in an enormous wheeling motion over the Gulf of Mexico and the Southern United States, first towards the West and then later towards the East as the moisture moves north of the 30th parallel.

The atmospheric water balance for a region, such as the State of Texas, can be computed by taking the vector field of atmospheric moisture flow, **q**, defined on the 2 degree grid, and interpolating the flow vectors onto points lying on the border of the State at regular intervals. The border of Texas is 7000 km long and points on the border at 100 km intervals seem to be appropriate. The lines joining these border points are also a set of vectors, **L**, which can be imagined to project vertically into the atmosphere and form a closed set of planes surrounding the atmosphere of the State, which are the control surface of an atmospheric control volume drawn over the State. The flow of water across any one of these planes can be found from the vector cross product:

Q = qxLy - qyLx

where the Easting and Northing components of the flow vector **q** = (qx, qy), and the border line vector **L** = (Lx, Ly) are cross multiplied. This is a rather interesting example of two different uses of the concept of a vector: firstly, in the sense of physical science for the atmospheric flow data, **q**; secondly, in the sense of GIS for the border line, **L**. By continuing this computation around the border, the net outflow of atmospheric moisture across the control surface can be determined, which is equal to difference between evaporation and precipitation on the land surface plus the change in atmospheric moisture content within the control volume. It turns out that the difference between atmospheric moisture inflow and outflow over Texas, Q, is usually less than 5% of the average atmospheric moisture flow over the State, so small errors in the flow computation can lead to substantial variations in the value of Q. The net evaporation from the land surface computed from this global data set may not be very accurate. The US National Meteorological Center is presently implementing a new mesoscale GCM over North America called the Eta model, using 40 km computational cells. The results of this model should permit more accurate atmospheric water balance computations to be made for regions like Texas.

Soil water is that water contained within the soil column, so the control volume for the water balance is a block of soil, and the computation consists in relating the change in soil moisture content to evaporation, precipitation, and outflow from the soil. At first sight, it would seem that the most appropriate spatial unit to use for a soil water balance would be a soil map unit, but these map units have very irregular shapes and a great range in size from one map unit to another. Climate data also play an important role and it appears that because the construction of spatially distributed climatic data usually involves interpolating climate onto a grid, that the grid cells used for that interpolation are also appropriate as "soil boxes" for constructing a soil water balance. Such climate cells are approximately 3' in size for the climate grids of the United States constructed by Daly et al. (1994), and for Africa constructed by Hutchinson et al. (1995). Willmott et al. (1985) constructed a global soil water balance model, and later Legates and Willmott (1990) developed a global climatology of monthly mean temperature and precipitation grids on 0.5 degree cells. Willmott uses the simple Thornthwaite soil water balance method, which requires only a single soil water parameter, the water holding capacity. Other soil water balance models, such as that in SWAT, have several soil layers and require more soil parameters. The recent emergence of a satellite derived net radiation balance of the earth (Darnell et al. 1992) provides net radiation estimates for the soil water balance, an important new data source.

The product of a soil water balance is a time history on a daily or a monthly basis of soil moisture content, evaporation and "water surplus" which is the water flowing from the soil to form surface runoff and groundwater recharge. Given the same input data, computation on a daily basis will always yield more water surplus than will computation on a monthly basis because daily precipitation is an episodic process, zero on most days, but when a precipitation event occurs, the soil moisture storage can be quickly filled up, thus producing a water surplus; if the same data are averaged over a monthly interval, it is as if the precipitation falls as a gentle mist, which may evaporate back to the atmosphere before the soil moisture capacity is filled. Interpolation of daily precipitation onto a grid is an uncertain undertaking because the spatial variation in daily precipitation is large. There is thus a challenge in constructing a GIS hydrology model for soil water balance in choosing the appropriate time interval for calculation.

There are two kinds of groundwater flow: *unconfined flow*, which occurs near the land surface and for which the water table or phreatic surface is the upper boundary of the saturated aquifer, and *confined flow* which occurs in deeper aquifers where the upper boundary on the flow is provided by an overlying aquitard, or hydrogeologic unit of low permeability. Unconfined flow is strongly influenced by surface hydrologic conditions, especially by seepage from streams crossing the region where the aquifer outcrops at the land surface. Confined flow is less influenced by surface conditions. Groundwater flow models usually use rectangular or triangular spatial units to represent elementary volumes of porous medium upon which the computations are to be performed. A GIS-based groundwater model can similarly be constructed using polygonal units provided the polygons are reasonably regularly shaped, but groundwater models constructed within GIS are not very computationally efficient.

In constructing a groundwater balance model, there are two computations to be performed: first, a water balance on each spatial unit in which all the inflows and outflows of the unit are used to determine the change in water storage and thus of the piezometric head within the unit; second, a flow computation between each pair of spatial units in which Darcy's law is used to determine the rate of groundwater flow as a function of the difference in head and the flow properties of the aquifer in the units. In a map-based groundwater modeling system, the first computation is done over all the polygons making up the aquifer, while the second is done over all the boundary lines of those polygons. Interaction between surface water in streams and underlying groundwater can be similarly determined by applying Darcy's law to the difference in piezometric head between the stream passing through an aquifer unit and the surrounding aquifer. All these computations need to be done on reasonably small units not more than say 20 km in cell size, because otherwise the head gradients in space become very small. Groundwater aquifers are usually quite confined in area and do not extend over the whole landscape, so unlike surface water flow which takes place everywhere, groundwater flow is more of a localized problem and a regional study needs to take into account each aquifer in the region individually, rather than considering groundwater flow to be a regional phenomenon.

Surface water is water in streams, lakes, wetlands and reservoirs. This water system is in some ways the most complex of all the phases of the hydrologic cycle, because it interacts with the other three phases, namely atmospheric water, soil water and groundwater, because the flow velocity is large compared to the velocity of groundwater flow, and because the flow environment is complicated, depending in part on the characteristics of the land surface and in part on the characteristics of the stream system. Fortunately, this is the area where GIS helps the most because of the detailed description of land surface features which can be presented in GIS. As described earlier, by making a suitable terrain analysis using DEM data, a conceptual model of the surface drainage system can be built up in which each watershed has one and only one stream draining it, and each watershed and stream pair can be assigned the same identification number. The watersheds so constructed are of two types: a source or head watershed in which the stream originates within the watershed, and an intermediate watershed where the stream flows both into and out of the watershed.

The stream network is manipulated so that each stream is represented by a single arc, and the arcs are flow ordered so that the from node is upstream and the to node is down stream. Each stream arc is enclosed within its associated watershed polygon. Watershed boundaries are delineated from each stream junction so at most a node can have two streams flowing into it and one flowing out of it. Three flow variables can be associated with each watershed: "From Flow", "To Flow", and "Polygon Flow". From Flow is that stream discharge at the from node; To Flow is the corresponding discharge at the to node; and Polygon Flow is that discharge which comes into the stream by drainage from the surrounding watershed. Polygon Flow is computed by applying a a unit hydrograph to the water surplus computed by the soil water balance model, and it may also include a component representing exchange of water between the stream and the underlying groundwater aquifer. This implies that the soil water surplus data may have to be spatially transferred from the soil water balance spatial units to the watershed units by using polygon overlay functions. The Polygon Flow can be divided by the length of the stream and added progressively to the discharge along the path of the stream, so in a time-averaged calculation, the discharge at a distance D from the upstream end of a stream of length L is given by:

Discharge = From Flow + (Polygon Flow / L) * D

In time-varying flow, the computation is more complex and stream routing methods such as the Muskingum method (Fread, 1993) are appropriate for computing the time distribution of the To Flow given the time distribution of the From Flow and the Polygon Flow. The time table structure described earlier in the paper is used to record the results of these calculations with a separate table being used for each of the three flow variables, a separate field for each watershed, and time on the vertical axis of the table.

In doing such computations, there is a choice between sequencing the computations "first in space and then in time" or "first in time and then in space", in other words, doing all the computations for a given watershed through time and then moving to the next watershed, or doing the computations for all watersheds for one time interval before moving to the next time interval. For regional hydrologic studies where the watersheds don't influence one another the "first in space then in time" method works best; for more localized hydraulic analyses where water level in the river influences groundwater levels, the simultaneous nature of the interactions makes the "first in time then in space" method necessary. The nature of the process interactions governs the computational sequence.

All of the preceding discussion is valid for a pristine landscape untouched by human activity. But reservoir construction, pumping of water from rivers and groundwater systems, and discharge of wastewater all have profound effects on surface and groundwater flow and quality. Modeling the effects of water utilization permits GIS hydrology models to be useful as a basis for planning decisions on water facilities. Consider a pumping station withdrawing water from a river. In a spatial sense, this can be represented by a point on a river with attributes describing the time pattern of withdrawals. It is useful to locate this point by its relative distance from the upstream end of the reach (D/L in the nomenclature defined previously). This "proportional aliasing" provides a way of locating objects on river reaches which is somewhat independent of the scale of the map used to define the river reach spatially. River reaches are always longer on maps of larger scale but the relative location of a point remains reasonably stable across scales. Once the withdrawal is located, it can be included in the discharge computation.

Including reservoirs in the computation is more complicated because a reservoir is a fairly complicated system by itself, and its outflow depends on the manner in which the reservoir control works are operated. A standard method of simulating water supply reservoirs is presented by Chow, Maidment and Mays(1988).

- The ready availability on Internet and CD-ROM of data describing the land surface, especially digital elevation data for land surface terrain, which has made it practical for the first time to delineate watersheds in a few minutes in an automated way, and to compute the hydrologic properties of those watersheds. The ability to build an integrated spatial hydrologic data base for a particular region has greatly improved.
- DEM data can be applied for small scale or large scale problems in a relatively scale-independent manner. The Digital Chart of the World, and the 30" DEMs derived from it by the USGS permit the extension to the world of grid-based hydrologic techniques.
- The interpolation of climatic data, especially mean monthly temperature and precipitation, has advanced to the point where its standardized data sets are very useful for regional hydrologic studies. In effect, the contour lines on the classic climatic atlas have been replaced by fine mesh digital grids of data.
- Comprehensive hydrologic simulation systems using GIS databases are now operational, and being used for analysis of basins of more than 1 million km2 in area, a task that would have been unthinkable a few years ago. Large regional hydrologic models are now within grasp.
- Some progress on integration of processes among different phases of the hydrologic cycle is being accomplished by formulating separate models for each phase and then using GIS spatial data handling capabilities to transfer results from one set of spatial model units to another.
- The advent of object-oriented GIS programming languages has broken the barrier to capturing time variation of spatial processes that was so great a limitation in earlier GIS applications to hydrology.

At the same time, a number of formidable challenges still remain before GIS achieves more of its potential in hydrology:

- The standardized data bases that are generally available are applicable to regional scale analysis of fairly large watersheds. Considerable database development is still needed within cities to support analysis at the scale of small urban watersheds. Probably, 80% of hydrologic modeling is done to solve problems in urban areas, so this limitation is critical.
- Various methods for creating GIS-based models of hydrologic processes are emerging but they have not yet been standardized to the point that they are being applied widely. There is a great need for dispersion of knowledge so that more people can use what is available. Internet will help with this. The emergence of more powerful PC-based GIS systems will also help.
- The integration of hydrologic processes, particularly integration of surface and groundwater flow, is not yet solved very well. Integration of processes across scales of space and time is not well understood. A map can be drawn at any scale, but it is unclear to what extent models can be applied at different scales.
- The impact of water utilization facilities, such as pumping stations and reservoirs, on flow through the landscape, is not well described in spatial hydrology models as yet.
- Subsurface representation of hydrogeologic properties of aquifers is embryonic. There are no standardized databases of hydrogeology, like the ones that exist for soils. A 3-D GIS system has not really emerged yet.
- Water quality modeling in rivers and lakes is sufficiently complex that it is still largely being done in traditional simulation models supported by GIS data. There is not yet much intrinsic water quality modeling within GIS.

But, looking back to the first conference four years ago, there is no question that significant progress has been made in spatial hydrology. Most of the challenges that I cited in my earlier paper are now covered in the accomplishments list of this paper; the set of challenges that I have identified here are largely new ones which hopefully will also largely be addressed during the next few years.

Chow, V.T., Maidment, D.R., and Mays, L.W. (1988) *Applied Hydrology*, McGraw-Hill, New York, pp. 530-537.

Daly, C., Neilson, R.P. and Phillips, D.L. (1994) A Statistical-Topographic Model for Mapping Climatological Precipitation over Mountainous Terrain, *Journal of Applied Meteorology, *33(2), 140-158.

Fread, D.L. (1993) Flow Routing, in *Handbook of Hydrology*, ed. by D.R. Maidment, McGraw-Hill, New York, pp. 10.5-10.13.

Hutchinson, M.F., Nix, H.A., McMahon, J.P. and Ord, K.D. (1995) *A Topographic and Climatic Data Base for Africa - Version 1*, CD-ROM available from Centre for Resource and Environmental Studies, Australian National University, Canberra ACT 0200, Australia.

Legates, D.R., and Willmott, C.J. (1990) Mean seasonal and spatial variability in gauge-corrected, global precipitation, *Int. J. Climatology*, 10(2), 111-127.

Maidment, D.R.(1993) GIS and Hydrologic Modeling, in *Environmental Modeling with GIS*, ed. by M.F. Goodchild, B.O. Parks, and L.T. Steyaert, Oxford University Press, New York, pp. 147-167.

Willmott, C.J., Rowe, C.M. and Mintz, Y. (1985) Climatology of the terrestrial seasonal water cycle, *Jour. of Climatology,* 5:589-606.

David R. Maidment

Dept of Civil Engineering

University of Texas

Austin Tx 78712

Telephone: (512) 471-0065

Fax: (512) 471-0072

Email: maidment@crwr.utexas.edu