Terrain-based catchment partitioning and runoff prediction using vector elevation data

WATERRESOURCES RESEARCH, VOL.27,NO.6, PAGES1177-1191, JUNE1991 Terrain-Based Catchment Partitioning andRunoffPrediction Using Vector Elevation Data IAN D. MOORE' Department of Agricultural Engineering, University of Minnesota, St. Paul RODGER B. GRAYSON Department of CivilandAgricultural Engineering, Centre for Environmental Applied Hydrology, University of Melbourne, Parkville, Victoria, Australia An automated methodof partitioning catchments intointerconnected elementsusinga "stream tube"approach andvectoror contour-based digitalelevation models is brieflydescribed. Withthis formof partitioning, hydrologic models canbe structured basedon thehydraulics of flowwithina catchment andtheeffectsof topography onrunoffproducing mechanisms andspatial!y distributed flow characteristics (suchasflowdepthandvelocity) canbedirectly,andrealistically, accounted forin the modeli.The method allowscomplex three-dimensional flowproblems to be reduced to a seriesof coupledone-dimensional problemsin areaswith complexterrain.Two simpleprocess-oriented hydrologicmodelsthat demonstrate the utilityof thisformof partitioning are presented. The first models subsurface flow-saturation overland flow and the second models Hortonian overland flow. Observedand predictedrunoffhydrographs and the dynamicexpansionand contractionof runoff sourceareason a smalllaboratorymicrocatchment are presented.Also shownare predictedrunoff hydrographs and surfaceflow velocitieson a smallrangelandcatchmentin the United States. INTRODUCTION Many hydrologicand water quality modelscrudelyrepresentthe three-dimensionalnature of natural landscapesand thereforecrudely representspatiallydistributedhydrologic processes.As transport modeling becomesincreasinglyimportant in hydrologic and environmental assessment, this becomesa limiting factor in the predictive power of these models. Not only do we now need to know the temporal variation in discharge at the catchment outlet, but we also needto be able to accurately predict the temporalvariation in flow depthsand flow velocitiesthroughoutthe catchment. For example,accuratemodelingof soil detachment,entrainment, transport and deposition throughout a catchment requiresa knowledgeof the temporaland spatialvariationin graphic attributes and terrain-based indices of a variety of hydrological,geomorphologicaland biological processesare discussedby Moore et al. [1991]. Readers are referred to this paper for a more detailed discussionof grid- and TIN-based DEMs. Vector DEMs are based on the most common form of elevation data storage, the topographicmap. Contour lines are represented in the vector DEM by a series of x, y coordinates (Cartesian coordinate system in the horizontal plane) for each isoline of elevation [Moore et al., 1988]. The spacingof the x, y points can be either regular or irregular with the latter requiringless data storagefor a given level of accuracy. Irregular spacing allows areas with relatively straightcontoursto be digitized at a larger spacingthan areas where the contour has more curvature. Topographic maps discharge,flow depthand flow velocity. The effectof topo- are prepared directly from aerial photographsor field surgraphicconvergenceand divergenceon flow characteristics veys so the information has undergone the minimum of in naturallandscapeshas a major impact on the valuesof manipulation,therefore minimizing errors. The major disadthesehydrologicand hydraulicvariables.Without the accu- vantagesof suchmethodsare that digitizingis relatively time ratepredictionof terrain characteristics,accurateprediction consuming and the data storage requirements are much of water quality is not possible. greater than for grid-based methods or TINs. The informaDigital elevation models (DEMs) provide the basic data tion density of vector DEMs is greater in areas with steep for characterizing the topographicattributesof landscapes. slopesthan in flatter areas, where the topographyis less well Thereare threeprincipalways of structuring a networkof defined. elevationdata: (1) triangulatedirregularnetworks,or TINs; In his critique of many physically based models, Beven (2) grid networks;and (3) vec,toror contourline based [1989, p. 170] argued that future developments"must take networks.The advantagesand disadvantages of thesethree accountof the need for a theory of the lumping of subgrid structures,togetherwith applicationsof terrain analysis scaleprocesses;for closer correspondencein scalebetween methods based on these structures for calculating topo- modelpredictions and measurements; (and)for close• correspondencebetween model equations and field process1Nowat the Centrefor Resourceand Environmental Studies, es.... "For these criteria to be achieved there is a need for Australian NationalUniversity,Canberra,AustralianCapitalTerri- a method of partitioning landscapesinto small areas where tory. the hydrologicprocessesand the soil, vegetation,and topoCopyright1991by the AmericanGeophysical Union. graphiccharacteristics can be considereduniformor at least can be characterizedby simple relationships.These must be •Paper number 9!WR00090. consistentwith the hydrologic processesoccurring in the 0043-1397/9 !/91WR-00090505.00 1178 MOORE AND GRAYSON: CATCHMENT PARTITIONING AND RUNOFF PREDICTION landscapeand provide a flamework for structuring both hydrologic models and geographicinformation systems, GISs, that are increasingly being interfaced with these models. With few exceptions(e.g., ANSWERS [Beasleyet al., 1980];AGNPS [Young et al., 1989]; SHE [Abbott et al., 1986]), hydrologic models are not structured in a way that allows them to easily accessand use data containedwithin a GIS. The ANSWERS, AGNPS and SHE models and most GISs use a square grid or cell network as their basic structure. However, the cell sizesused in a GIS and hydro- logic model are often inconsistent,and meaningfulaggregation and disaggregationof data from one cell size to another is not a trivial exercise. Mark [1978, p. 24] notedthat grid structuresfor spatially partitioningcatchmentsare not appropriatefor many geomorphologicaland hydrologicalapplications.He statedthat "the chief sourceof this structureshouldbe the phenomena in question, and not the problems, data, or machine considerations,as is oftenthe case." For hydrologicmodelingand analysisin complexterrain, the processin questionis the hydraulicsof fluid flow in the landscape.Vector DEMs are capableof providinga structurebased on the physics of water movement on the land surface without requiring a priori knowledgeof the terrain [Moore and Grayson, 1989; Moore et al., 1991].Grid structuresdo not meetthis requirement. They typically restrict flow from a given node to only one of eight possibledirections,and computedflow paths take on a zigzagshape.Recently, Vieux et al. [1988]useda TIN to providea frameworkfor solvingthe kinematicforms • Defined Catchm•_nt boundary Contour -- -- -- lines Streamlines Fig. 1. Hypothetical catchmentshowingthe elementsformed of the two-dimensional overlandflow equationsusingthe by the intersectionof adjacentstreamlinesand equipotentiallines finite elementmethod.In this applicationthe quadrilateral (contourlines).The shadedregionis a singlestreamtubeformedby and triangularelementswere alignedso that the principal a pair of adjacentstreamlinesand consistsof several serially slope direction was parallel to one of the sides of the connectedelements(adaptedfrom Ohstad [1973]). element.This satisfiesMark's criteriafor a suitablehydrologicallybasedstructure,but it requiresa prioriknowledge of the topographic attributesof the landscapeto suitably then assumedto have uniformproperties.We then demonalignthe elements.Vieuxet al. [1988]developed theTIN by stratehow this form of partitioningcan be usedto structure fitting it to a vector DEM. dynamichydrologicmodelsfor determiningthe temporal The first vectoror contour-based methodof partitioning and spatial variation in hydrologicresponseto climatic catchmentswas proposedby Ohstadand Brakensiek[1968] inputs. Two simple kinematic models were developed who termedit a "streampath" or "streamtube" analogy. around this topographicstructure. They simulate the dySpeight[1974,1980]subsequently usedthe conceptto clas- namic and spatiallydistributedhydrologicresponseof a sifylandforms, especially specific catchment areas(upslope catchment where (1) subsurface flow and saturation overland contributingarea per unit contourwidth) and specificdis- flow runoff mechanismsdominate, and (2) Hortonian overpersal areas (downslopedispersalarea per unit contour !and flow dominates. width). Specificdispersalareascannotbe computedwith current grid-basedmethodsbut can be easily calculated CATCHMENT PARTITIONING USING TAPES-C usingthe contour-based approach.Grid-basedmethodsof structuring andanalyzing terraindatatreatrunofffroma grid In using Onstad and Brakensiek [1968] "stream tube" cell(twodimensional) asa pointsource(zerodimension) and concept,contourlinesare assumedto be equipotentiallines, projectitsmovement downslope by a line(onedimensional). and pairs of orthogonalsto the equipotentiallines (streamHowever, the vector-based structureusingthe "stream lines)form the "streamtubes" (Figure 1). Cayley[1859,p. tube" methodof analysisis morephysicallyrealisticas it 264] called thesestreamlines"slope lines (lignesde la plus treatsthe water as originatingfrom a segment(one dimen- grande pente)" and "orthogonal trajectories," and desional)of an equipotentialline and projectsonto a two- scribedthe generalrelationshipbetween contour lines and dimensional surfacedownslope. It shouldalsobe notedthat slope lines at summits, immits (depressions),ridges and oncethe partitioning of a catchment is completeusingthis saddlepointsin the landscape.Adjacent contour lines and "streamtube"technique, thedatastorage requirements are streamlinesdefine irregularly shaped elements(Figure 1). no more onerous than for any other form of DEM because Surfacerunoffentersan elementorthogonal to the upslope only the element attributes must be stored. contourline and exits orthogonalto the downslopecontour The first part of thispaperbrieflydescribesan automated line, with the adjacentstreamlines beingno-flowboundaries. computer-based methodof partitioningcatchments into ele- Flow from one elementcan then be successivelyroutedtO mentsusingthe "stream tube" concept.Each elementis downslopeelementsformed by the same streamtube. With MOOREANDGRAYSON' CATCHMENT PARTITIONING ANDRUNOFFPREDICTION thisform of partitioning a catchment, one-dimensionalflow 1179 I canbe assumedwithin eachelement,allowingwatermove- mentin a complexthree-dimensional catchment to be represented by a seriesof coupledone-dimensional equations. Theequationscan be solved usinga one-dimensionalfinite differenceor finite element scheme. TAPES-C, TopographicAnalysis Programsfor the Envio ronmental Sciences-Contour, is a setof computer programs thatautomatically partitiona catchment intointerconnected, irregularly shapedelementsusingthis "streamtube" approach. The methodis a development of the contour-based HP3 terrainanalysismethodsdescribedby Moore et al. [1988] and O'Loughlin [1986]. Application of these earlier terrain analysis methodswas restrictedto steadystatehydrologic modelingor quasi-dynamicmodeling(as used by Moore et al. [1986])becausethe contour lines and computedstreamlines did not form a series of interconnected elements that would allow flow to be routed from one element to another. TAPES consists of four programs: (1) DIGITIZ, which permitscontourlinesto be line-digitizedto createa vector DEM; (2) PREPROC, which is a preprocessing programthat transforms the vector DEM into a north-south oriented coordinate system; (3) TAPES-C, which subdivides a catch- 0 250 500 I • I SCALE $ mm Fig. 2. Topographic map of the 2.0-m2 sandbed microcatchment (10 mm contour interval) showingits subdivisioninto elements using TAPES-C and a saddle point (S) and a poorly defined interfluve (I). mentinto elementsusing the "stream tube" approachand calculatesa variety of topographic attributes for each element;and (4) TAPES-G (TAPES-grid) which is a grid-based methodof terrain analysis that will not be discussedhere. that only the x, y coordinates of points defining two contour The main functions of programs DIGITIZ and PREPROC have been described previously by Moore et al. [1988]. lines need to be stored in computer memory at any one time However, they have been modified to accept a larger data rather than the x, y coordinates of all contour lines as in the earlier modelsof Moore et al. [1988] and O'Loughlin [ 1986]. set, accommodate a variable contour interval so that increasedtopographic detail can be input where necessary, Furthermore, only four points are needed to define each andpermit the catchment boundary and the x, y coordinates element. This reducesthe computer memory requirement of of points of interest on the catchment to be input directly the program and allows TAPES-C to be run on small during digitizing. Furthermore, the Universal Transverse personal computers with limited memory capacities. The Mercator Cartesian coordinate system is now used as the applicationof TAPES-C to a small catchmentis illustrated in basefor all calculationsrather than rescalingthe coordinates Figure 2. The spacing between adjacent streamlines, X, and hence of the vector DEM to lie within the bounds of a 1000 x 1000 arbitrary unit coordinate system. Program PREPROC can the maximum size of each element or the width of a "stream alsogenerate a vector DEM in the form of sets of x, y, z tube," is controlled by an input parameter AINT. If 1.3AINT coordinatesthat can be subsequentlyused with other inter- ß< •. < 2AINT a new streamline start point is inserted midway (X/2) between the existing streamlines on the downstream polationprograms to generate a grid DEM. Moore et al.'s [1988] approximate method of estimating contour line, as indicated by A in Figure 3. If 2^Isz < X < flow trajectories is used in TAPES-C to determine the n ^INT, where n is an integer greater than 2, then n - 1 equally spaced [X/(n - 1)] streamline start points are streamlines or "stream tubes." The streamlines between contourlines are approximatedby short straight-lineseg- ments,which simplifiesthe computer program and minimizesthe computationaltime requiredto subdividea catchment. Streamlinesare computed using two criteria: (1) minimum distance between adjacent contour lines [O'Loughlin,1986]; and (2) orthogonalsto the downslope contour lines. These two criteria are used in an attempt to overcomethe error causedby usingstraight-linesegmentsto define streamlines. The first criterion is used in ridge areas andthe secondcriterion is used in valley areas.Application of the two criteria is determinedby the curvatureof the 25 m 25 m 20m'" ••.•••y,• ",••\ 20m contour lines(i.e., plancurvaturein radianspermeter).The 15 rn '" •.•.• 'V' 'NN• • !5 m contourcurvatureis computedby numericaldifferentiation 10m •10m of nodes(of x, y coordinates) definingthe contourline. TAPES-C performsthe partitioningof the catchment Fig. 3. Subdivision of a series of hypothetical contour line beginning at the contourline of lowestelevationandending segmentsusingthe minimum distanceand orthogonalcriteria and at the highestcontourline, successively determiningthe showingtrajectorystartpoints(A and B) and a trajectoryendpoint elements for eachadjacentpairof contourlines.Thismeans (C). 1180 MOOREAND GRAYSON:CATCHMENTPARTITIONINGAND RUNOFF PREDICTION insertedbetweenthe existingstreamlines,as indicatedby B in Figure3. It canbe seenfrom Figure2 that an elementcan havemorethanoneupslopeconnecting element,but thatan proposed by Band[1986].A streamchannelpasses through streamlines at a contour line when ;t < AINT/2, as illustrated ture of the contourlines.The streamchannelor waterway an element if the total upslope contributing area of the elementis greater than a critical value. The streamchannel elementonlyflowsintoonedownslope element.In an earlier is assumedto intersect the upslope and downslopecontour version of TAPES-C, we included the option of stopping lines bounding an element at the points of maximum curva. by C in Figure3. Thisproducedan elementnetworkthat was shownin Figure 4 was computedby TAPES-C usinga uncluttered and elegant in appearance.However, when critical upslopecontributingarea of 4.2 ha.. usingtheresulting elementnetworkfor runoffprediction,we couldnot developa reliablerulefor distributing the outflow MODELING RUNOFF PRODUCING MECHANISMS USING to multipledownslope elements.Consequently, this option TAPES-C was eliminated. The ability to vary the contour interval in TAPES-C allowsthe userto input a densersamplingof contourlinesin regions of flat topography. However, complete contour lines, extendingfrom one catchmentboundaryto another, mustbe inputin theseregionsratherthan shortcontourline segments.The previousmodelsof Moore et al. [1988] and O'Loughlin [1986]had difficultyproducingrealisticstreamlines in flat regionswhen the contour lines were far apart becauseof the assumptionof straight-linesegmentsto define Ohstad and Brakensiek [1968], Kozak [1968], and Onstad [1973] used the "stream tube" method to examine the distributedrunoff behavior of small catchmeres.However, in all four studiesthe catchmentswere partitioned by hand. The studies of Onstad and Brakensiek [1968] and Kozak [1968] were for small catchmentswith very simpletopographies, while that of Ohstad [1973] was for the simple hypothetical symmetric catchment shown in Figure 1. Recently, Tisdale et al. [1986] used this technique togetherwith streamlines between adjacent contour lines. Varying the an implicit one-dimensional kinematic wave equation to contour interval overcomesthe requirementof selectingthe examine distributed overland flow and achieved good accucontour interval of the input topographic data to fit the racy when comparedto steady state solutionsfor flow depth and discharge over a hypothetical cone-shaped catchment. flattest part of the terrain. Current topographicmaps generally use a constant con- They found the method to be computationally faster than tour interval and do not define certain types of topographic two-dimensional finite element modeling methodologies,besingularitiessuch as pass or saddlepoints [Warntz, 1975]. cause it solves simpler model equations, while retainingan Furthermore, becauseof the largecontourinterval, typically equivalent physical realism [Moore et al., 1991]. We will now demonstrate how TAPES-C can be used as a 20 ft (6 m) for the 1:24,000 scale United States Geological Survey map seriesand 10 m or 20 m for the 1:25,000scale framework for distributed hydrologic modeling that accounts Australian map series,ridgesand interfluve areas, particu- for the effect of three-dimensional terrain on storm runoff larly at or near the catchment boundary, are often poorly generation.We present two simplified,but physicallyrealdefined. For example, in the Piedmont region of the eastern istic, process models: one based on the mechanismsof United States, topographicsingularitiessuchas peaks often saturation overland flow and subsurface flow, and the other do not occur, or at least cannot be easily defined. Rather, based on the mechanism of HortonJan overland flow. Both level interfiuves, ridges or tablelands occur. In the region models use the kinematic approximationsof the governing where such topographic features appear, streamlines are flow equations.A terrain-basedmodel that integratesall difficult to estimate with the limited information appearing three storm runoff generatingprocessesis the subjectof a on the topographicmaps.TAPES-C estimatesthe streamline separatestudy [Grayson, 1990]. in such cases by projecting an orthogonal from the downslope contour line upslope until it intersectsthe catchment Basic Stream Tube Structurefor Hydrologic Modeling boundary. The downslope contour line, the two adjacent The TAPES-C method of partitioning natural landscapes streamlines and the catchment boundary then define an element. Figure 2 demonstratesthe application of this ap- into a series of interconnected elements allows the twoproximation in the region of a saddle (S) and a poorly dimensionalhydraulic equations describingwater movement to be simplifiedto a seriesof defined interfluve (!). This method of representationalso on and in suchlandscapes allows TAPES-C to be easily applied to small-scale plot coupled one-dimensionalequations. Figure 5a presentsa segmentof a catchmentand showshow it mightbe subdistudieswith imposedboundaries. For each element, TAPES-C computes the following vided into elements using TAPES-C. Figure 5b showsan networkto whichthe topographic attributes: (1) element area; (2) total upslope equivalentcoupledone-dimensional contributingarea; (3) connectivityof upslopeand downslope one-dimensionalfinite difference or finite element forms of elements; (4)x, y, z coordinatesof the element centroid; (5) the governinghydraulicequationscan be applied.Nodal to the midpointson the x, y, z coordinates of the midpoint on the downslope pointson the networkcorrespond contour bounding the element; (6) the average slope of the upslopeand downslope contourlines boundingeacheleelement orthogonal to the contour; (7) the widths of the ment.Thegoverning flowequations aresolvedat eachnodal element on the upslopeand downslopecontourlines bound- point in the network. ing the element; (8) the flow distanceacrossthe element; and This methodof solvingthe flow problemand structuring (9) aspect or azimuth of the element. The computation of the networkrequiresestimatesof the plan area of each these attributes has been previously describedby Moore et element, Ae, theflowdistance alonga streamline through al. [ 19881. Streamchannelnetworksare definedby TAPES-C usinga method based on the upslope contributing area criteria the element,ASe, local slope,/3, the widthsof the flow domain(i.e., of the elements),w, at eachnodeand the connectivity of upslope anddownslope elements asbasic MOOREANDGRAYSON: CATCHMENT PARTITIONING ANDRUNOFFPREDICTION tl I1! III 1181 11t HP2 111 N Metres 5O I00 '/" • , "• 200 , 400 Feet Wolerwoy Fig. 4. Subdivisionof the Chickasha,Oklahoma,R-5 catchmentusingTAPES-C. Shownis the digitizedvector DEM (the 0.6-m intervalcontourlines),digitizedhighpoints(HP1, HP2), catchmentboundary,locationsA, B, and C at which runoff hydrographsare predictedand the predictedlocationof the streamchannel. inputparametersto the hydraulicequations.Each element hasa uniquenumberand can be connectedto any numberof upslopeelements, but only one downslopeelement. The elementof highestelevationis numbered1 and the numberingproceedsalongeach contourline in a consistentleft-to- right or right-to-leftdirection and endswith the last element on the contour line of lowest elevation. When two upslope elements flow into one, the flows are combined and used as inputto the upslopeface of the singleelement,as shownfor the lowest element in Figure 5a. TAPES-C automatically providesestimatesof these attributesfor each element as it carries out the catchment partitioning. Adjacent elements in depressionalareas with convergent flow lines can have widely varying upslope contributing areas so that the distribution of flow cross-sectional areas and dischargesacrossthe upslope face of an element with a large number of tributary (i.e., upslope) elements can be highly variable. In the routing algorithms, the dischargeand flow cross-sectional Contour lines a) Streamlines b) area into an element is used in the calculation of outflow discharge and cross-sectional area. Hence the way the elements are connected for flow routing is important. Three approacheswere examined for determining the inflow discharge and cross-sectional area to an element with multiple tributary elements: (!) sum the tributary outflow cross-sectional areas to provide the inflow cross-sectionalarea to the element; (2) sum the tributary dischargesto the element and calculate an equivalent flow Fig. 5. (a) A hypothetical hillslopesegment 'subdivided into cross-sectional area based on the propertiesof the element; elements usingTAPES-Cand(b) its equivalent one-dimensional network. and (3) assumeflow is channelizedand that only the upslope 1182 MOOREAND GRAYSON:CATCHMENTPARTITIONINGAND RUNOFF PREDICTION tributaryelementwith the dominantdischargecontributes of catchmentoutflowrather than distributedrunoff [e.g., Lee directlyto inflowandthat the otherelementsbecomelateral and Dellcur, 1976;Bevenand Kirkby, 1979].Also, many inflow to the channel as it passesthrough the element. Methods 1 and 3 producedsimilar resultsbut method 2 models only simulate the response of nonconvergent hillslopes, thus only partially representing the influence of produceda morerapid hydrographrise and fall. Method2 topographyon subsurface stormflow [Freeze,1972;Beven, producesa smallerinflow area than the other methods 1981,1982a,b; SloanandMoore, 1984].Convergent topogbecauseof the nonlinearityin the cross-sectional area versus raphy results in the concentration of subsurface flow in dischargerelationship.Thereforethe total inflow area re- certain areas of a catchment, which increase the likelihood quiredto producethe sumof the discharges is lessthanthe of surface saturation as the capacity of the soil profile to averageof the tributary inflow areas. Mass continuitywas transmit the water is exceeded [O'Loughlin, 1981; Beven maintained with all three methods for surface flow routing but was violated for combined surface and subsurface flow with method 2. Method 2 producesseriousoscillationsand numericalinstabilitiesin the simulateddischargesbecauseof the form of averagingpreviouslydescribedand becauseof the large changesin the cross-sectionalarea versus discharge relationship that occur at saturation when some tributaryelementsmay be eithersaturatedor unsaturated.In the models describedbelow, method 1 was used unless the element was a channel element (i.e., had a defined channel passingthroughit), in which casemethod 3 was used. and Kirkby, 1979]. In this way the influence of subsurface storm flow on the storm runoff can be greatly enhanced. Smith and Hebbert [1983], Takasao and Shiiba [1988], and Sunada and Hong [1988] developed kinematic surfacesubsurface flow models that account for flow convergence and divergenceon simple topographiesthat can be approximated by planar, conical and partial conical surfaces.The TAPES-C analysisallows the modelingof irregular topographies that cannot be representedby these idealized, smooth surfaces in a similar way. Hortonian overland flow. The Hortonian mechanism of runoff generation is most important in semiarid and desert Storm Runoff-Producing Mechanisms areas [Freeze, 1974; Sivapalan et al., 1987] and on agriculThree principal mechanismsproduce storm runoff: (1) tural land when surface sealing occurs. It has formed the saturation overland flow that occurs when rising water tables basis for the hydrologic component of most soil erosion intersect the soil surface, generatingexfiltration; (2) Horto- models [e.g., Beasley et al., 1980; Knisel, 1980]. Accurate nian overland flow that occurs when the rainfall intensity prediction of the Hortonian response of a catchment relies exceeds the infiltration rate of the soil; and (3) subsurface on a thorough understanding and effective representationof flow in which water flows laterally through a highly conduc- the infiltration processwhose literature is extensive and has tive soil profile to the stream [Dunne, 1983; Smith and been reviewed by Skaggs and Khaleel [ 1982]. However, two Hebbert, 1983; Freeze, 1972; Hewlett and Hibbert, 1967; major problemsexist: (1) the derivation of field valuesfor the Horton, 1933]. These mechanismsare part of a continuumof relevant parameters [Knapp, 1978], and (2) characterization processesand may operate singularly, but more commonly of the effect of spatial variability in infiltration on runoff in combination [Freeze, 1972]. In the cases of saturation [Sharma et al., 1980; Lane et aI., 1988; Loague, 1988]. The overland flow and Hortonian overland flow, precipitation desire to use distributed parameter models has exacerbated falling directly on the saturated zone at the soil surface the difficulty of parameter estimation by requiting a large produces surface runoff or overland flow. These saturated number of values rather than relying on values that characareas may occupy only a portion of a catchmentand may terize large areas [Quimpo, 1984]. vary in size dependingupon soil propertiessuchas saturated Despite these difficulties, several models of infiltrating hydraulic conductivity, the moisture release curve, organic overlandflow have been developed,includingthe numerical content, depth to restricting layers, antecedent soil water models of Smith and Woolhiser [1971], Smith and Hebbert content, and topography. Elsewhere, all the precipitation [1983], and Abbott et al. [1986], and the analytical solution infiltrates. Runoff generated in these zones is referred to as by Woodsand Ibbitt [1988]. The model presentedhereis partial-area or variable source area runoff [Betson, 1964; similar in many respects to these models and includes Hewlett and Nutter, 1970; Dunne and Black, 1970; Dunne, algorithms for simulatingHortonian overlandflow, as either 1983]. broad sheet flow or flow in microchannels, and flow in Saturation overlandflow and subsurfaceflow. Subsurface storm flow is generally considered to occur as lateral movement of water in the upper soil layers. Van de Griend and Engman [1985] reviewed the reasons for this and reported the influence of hard pans and impeding layers channels and streams. The element structure of the model, derived from TAPES-C, allows each element to have different infiltration characteristics and parameters. However, asnotedabove,to do thisis usuallynot practicalbecause of In mostof our applications of [Weyman,1973],ploughpansin agriculturalareas,soilpipes the largedatarequirements. [Tanaka, 1982; Tanaka et al., 1988;Jones, 1981]and macro- the model we measurethe infiltrationparametersfor each pores [Fletcher, 1952; Mosley, 1979] or highly permeable soiltype on a catchmentandmakethe crudeapproximation upper soil horizons [Whipkey, 1969]. When subsurfaceflow that theseparametersdo not vary within eachsoiltype.We converges, the capacity of the soil to transmit the flow is then numericallyoverlay a soil survey map over the element exceeded and saturated areas are formed. These saturated mapof the catchmentand assigna soiltype to eachelement areas are effectively impermeableso in additionto exfiltra- basedon the soiltype at the centtoldof the element.This approachassumes thatthe variationin infiltrationproperties tion, all rain falling on them becomesrunoff. Most attemptsto model the processesof subsurfaceflow withina soiltypeis lessthanthatbetweensoiltypesandthat and saturationoverlandflow have beenrelatively simplistic the within-elementcharacteristicsare uniform. Neither of [Boughton, 1987] and have been directed toward estimation these assumptionsis necessarilytrue [Beven, 1989]. MOOREAND GRAYSON: CATCHMENT PARTITIONING ANDRUNOFFPREDICTION 1183 model it is assumedthat the entire soil profile first wets up to i fieldcapacity.Followingfrom the resultsof Andersonand Burt's [1977]study,the capillaryfringe is not differentiated from the saturated zone. The resistanceequationsor kinematic forms of the momentumequationsfor subsurfaceflow only, and for combined subsurface flow and saturation overland flow, can be expressedin terms of' the Darcy and Manning equations. They are, respectively, A Element Boundaries Q=K--sin/3=roHKsin/3 ¾ Fig. 6. Schematicrepresentationof the flow systembeingmodeled by the simplified subsurfaceflow-saturation overland flow 0--<A<•o•,D or 0 <- H < D (3a) A > •o7D (3b) model. Q = toDK sin/3 + aA•m where K is the effectivehydraulicconductivityof the soil,/g Saturation Overland Flow and Subsurface is thelocalslope(in degrees), a = n-lto-2/3tanl/2/3and Flow Model The subsurface flow and saturation overland flow model m = • forbroad sheet overland flow,a = n- ••r2/3 tan1/2/3 andrn = ½for channe!ized flow,n is Manning's roughness presentedhere is based on the assumptionthat infiltrated coefficient,and •ris the coefficientin the assumedrelationwater flows downslope in a saturated layer overlying an ship between hydraulic radius, R, and cross-sectionalarea impermeablebase. If the subsurfaceflow rate exceedsthe of a channel, R = •Ac•/2 [MooreandFoster,1990;Moore capacityof the soilprofile to transmitthe water thensurface and Burch, 1986]. saturation occurs and rain falling on the saturated area High hydraulicconductivitiesare indicatedby the small becomes direct runoff. It is assumed that the infiltration responsetimes of many catchments to subsurface flow, capacity of the soil is always in excess of the rainfall intensity,so that Hortonian overland flow is ignored.Hence surface runoff can only be generated from rain falling on areas that are saturated and from exfi!tration of the subsur- face flow. The model was kept as simple as possible,yet is physicallyrealistic. We have adapted the combined surface-subsurface kinematicmodelingapproachdescribedby Takasaoand Shiiba [1988],that they attributed to Takasao et al. [1976], which doesnot require the explicit calculationof exfiltration.The one-dimensionalform of the continuity equation that is appliedto each element can be written as OA •+ Ot OQ •-- Os iAe ASe (1) where A is the "apparent" cross-sectionalflow area [Takasaoand Shiiba, 1988], Q is the discharge,i is the rainfallintensity,A e is the planareaof the element,andAse is the flow distancealonga streamlinethroughthe element. The "apparent" cross-sectional flow areasfor the two cases of subsurfaceflow only, and combinedsubsurfaceand saturationoverland flow are, respectively, A=•oyH 0-<A<•o3•D or0<-H<D (2a) A = •oTD+ A c A > •oyD (2b) where •o is the width of the element orthogonal to the streamlines, y istheeffective porosity (totalporosity - field particularlyforestedcatchments.This may be due to water flowing through interconnected macropores created by roots, other organic matter and cracks between aggregates and rocks rather than through the soil matrix [e.g., Mosley, 1979; Sloan and Moore, 1984]. Hence the hydraulic conductivity used to characterize the subsurface flow may be different from that measured from small soil cores. For field applicationsof the model this parametershouldbe considered as an "effective" hydraulic conductivity of the soil profile rather than the measuredvalue of hydraulic conductivity for the soil matrix. Brakensiek's [!967] four-point infinite difference solution of both the kinematic overland flow and channel flow equations is used to route the subsurface and surface flow between elements. This method of solution was used prima- rily becauseit is unconditionallystablebut also becausethe dispersionof the numerical schemeovercomesthe problems of kinematic shocks that can occur on cascading planes when solving the kinematic equations using the method of characteristics. However, Kibler and Woolhiser [1970] have shownthat numericaldispersioncausesan under prediction and delay in the magnitudeand timing of the peak runoff with this method and is a recognized limitation. The finite difference form of the continuity equation can be expressed as At At At +2iAe A;' e A4+2Q4 •se=Al-A2+A3+2Q2 ASe (4) capacitysoil water content),D is the thicknessof the where the Q and the A variables are the discharges and hydrologically active soil profile(abovean impermeable"apparent" flow cross-sectionalareas, respectively, At is layer),andH is the depthof flow abovethe impermeable the time interval, subscripts 1 and 3 refer to time t at layer,as shownin Figure6. In (2b), A c = to(H - D) for positionss and s + Ase, respectively, and subscripts2 and broadsheetflow,whereH - D is the overlandflowdepth 4 refer to time t -e At at positions s and s + Z•se, (H >D), andAc is thechannel cross-sectional areawhen respectively.At any time and for any element the terms on the flow concentratesin fills or defined channels.In the the right-handsideof (4) are known quantities,and (4) can be 1184 MOORE AND GR•YSON: CATCHMENT PARTITIONING AND RUNOFF PREDICTION solvedfor A 4 usingthe Newton RaphsonmethodsinceQ4 can be expressedas a function of A 4. Boundaryconditions. The boundaryconditionsfor the model were adopted to approximatethe conditionsthat would occur on a hillslopewith an impermeablelayer at the base of a relatively thin hydrologicallyactive zone. In the 50 '"----' --- 40 er- 20 measured examplegiven below we assumethat the slope of this impermeableinterfacecan be approximatedby the land surface and that the thickness of the hydrologically active zone is constant. Neither of these assumptionsis necessary if sufficient data are available to characterize their spatial 10 properties.However, such data are rarely available.Similarly, a seepageboundary could be defined instead of a no-flow boundary at the base of the hydrologically active zone if the site data permittedit and the continuityequation (1) were modified accordingly. At the upslopelimit of a streamtube, which may be a high point or a catchment divide (i.e., the first node in the network where s = 0), the model assumesa zero flow depth 0 0 20 40 60 80 1O0 120 140 Time (min) Fig. 7. Predicted and observed total runoff hydrographsfrom the sandbed catchment. condition: H(s, t) = 0 Q(s, t) = 0 at s = 0 (5) where s is the distancedown a streamtube from the point of highestelevation. Although the kinematic wave equationsdo not require a lower boundary condition, a seepage face occurs at the catchment outlet (i.e., at s = S). The rainfall intensity is assumedto be temporally variable but spatially uniform. Model simulation. The modelwas evaluatedon a 2.0-m2 sandbedmicrocatchment [Grayson, 1990]. The sandbedsurface represented topographic convergence and divergence and produced saturatedsubsurfaceflow and variable source area runoff, both of which are evident in the real world. The surface configurationcontained steep, short slopesand the sand had a high hydraulic conductivity, providing an extreme test of the model. A sculptured impermeable surface was first formed to representa complex topographywith convergingand diverging segments,and then a uniform thickness (65 --- 6 mm) of wet fine sand was packed on top of this surface. A topographic map of the catchment is presentedin Figure 2. This map was computed using stereo photogrammetry with a standard error in elevation of less than 2 mm. Rainfall was applied to the microcatchment by an oscillating nozzle rainfall simulator at a rate of 55 --- 5 mm h -•. The kinetic energy of the simulatedrainfall was low enoughnot to cause soil splash on impact. Surface and subsurfacerunoff were independentlymeasuredand the saturatedzoneswere photographed during the runoff event and later digitized and rectified. The saturated hydraulic conductivity of the re- packedsandbedwas 1530___ 230 mm h-• (five measurements)andthe totalporositywas0.38 _-+0.01 m3 m-3 (14 of rainfall. The first lasted 31 min and was followed 35 min later by a 21-min burst. Figure 7 shows that the model performs well on the rising limbs of the hydrographsbut not as well on the recession limbs. The division between surface and subsurface flow is presented in Figure 8, and it is clear that surface runoff rates are consistently overpredicted and subsurface runoff rates are underpredicted. This is partly due to the difficulty of matching the model boundary conditions to the prototype. The downslope boundary for subsurface flow in the prototype was a seepage face in which the potential gradientswere steeper than assumedin the model so that subsurface flow was underestimated. The model assumesthat the equipotential lines can be represented by surface contour lines. In an attempt to overcome this difficulty the results presented in Figure 8 assume that any surfacerunoff predicted for the elements in the region of the lower boundary, except at the outlet for surface runoff, would actually be subsurface runoff in the prototype. The dynamic nature of the predicted and observed areal extent of the zones of surface saturation at three different times during the first rainfall period is presentedin Figure 9. The general trends are modeled well but the overall size of the saturated zones are underestimated. There are two 50[ - / f =• ! I ! ! i I • measured subsurface simulated surface ated subsurface measurements). Measurements of the soil water retention curves and tests performed on repacked sand columns indicated that hysteresis occurs in the sand and that a distinct capillary fringe exists after initial wetting. The thickness of the capillary fringe was about 50 mm. The maximumdrainableporositywas approximately 0.20 m3 m-3 andthe drainable/tillable porosity in theregionof the capillaryfringefollowinginitialwettingwasabout0.08 m3 m-3 [Grayson,1990]. The simulatedand observedtotal runoff hydrographsare shown in Figure 7. The rainfall event consistedof two bursts õ 20 lO 0! 0 ' ' 20 ,I •--:."",---"'"'---w--!. , t.•'--•• 40 60 80 1 O0 120 _ , 140 Time (min) Fig. 8. Predictedand observedsurfaceand subsurface runoff hydrographsfrom the sandbedcatchment. MOOREANDGRAYSON: CATCHMENT PARTITIONING ANDRUNOFFPREDICTION • t LL • -!-IP1t / 1185 300 290 •' N--'--Hydraulic Head • E,•280 270 'HP3 uJ260 250• 350 45O ' ' 55O Positionalongtransect (mm) ! ' 650 B' Fig. 10. Predicted hydraulic head at 23.5 min and soil surface elevation along the transect B-B' shown in Figure 9c. ing areas while those immediately adjacent may have low upslopecontributingareas. Figure 10 presentsthe predicted hydraulichead acrossa sectionB-B' of Figure 9, orthogonal to the principalflow direction down the valley at the time of peak flow. It showsthat significantsubsurfacewater ridging is predictedwhere water concentratesin the valleys creating large hydraulic gradientsaway from the valley. This tends to dispersewater away from the drainagelines, increasingthe size of the saturated zones, as observed in the prototype. However, the model assumesthat the land surface gradient representsthe subsurface hydraulic gradient and so can not representthe broadeningof the saturatedzone caused by the ridgingeffect. The high saturatedhydraulic conductivitiesin the prototype create a worst case scenario for the occurrence of this phenomenon. In most natural landscapes, hydraulicconductivitiesand slopesare smaller and flow path lengths are longer, thus damping out this effect. However, groundwater ridging has been observed in some natural HP3 b \ ,, catchments. \ Hortonian • • • LL Overland Flow Model In the Hortonian overland flow model the infiltration rate, infiltration volume and rainfall excessare computed for each element using a relationship based on an equivalence between the Green and Ampt and the Horton infiltration capacity equations [Morel-Seytoux, 1988a, b]. Spatially and temporally variable rainfall intensity and infiltration is simulated, but the rainfall intensity in any time increment, At, is assumedconstant. The rainfall intensity and infiltra- 'HP3 tion behavior is assumed to be uniform on each element but varying between elements. The time from the onset of rainfallto surfaceponding,t•,, is givenby K(Os -- Oi)H c t•,= i(i- K) (6) -. c Fig. 9. Predicted(stippledarea) and observed(boundedby dashedline)surfacesaturatedzonesat threedifferenttimesduringa rainfallevent:(a) 2.5 rain, (b) 4.25min, and(c) 23.5min. where •cis the resaturatedhydraulic conductivity, Osis the field saturated water content, 0i is the initial soil water content, Hc is the effective capillary drive, and i is the rainfall intensity [Mein and Larson, 1973; Morel-Seytoux, 1988a].Aftersurfaceponding occurs,i.e., for t > t•,, the infiltration capacity is calculated by 1 reasonsfor thesedifferences:(1) the effectof the capillary f= K+•-.•(ip-g)(e kat-1)e -k(t-tp) (7) fringewasignoredin thisrelativelysimplemodel,and(2) subsurface waterridgingoccurredin the valleys.The ele- wheref is the average infiltration rate between times t - At mentsin thedepression areashavelargeupslope contribut- andt, i•, is therainfallintensitythatinducedponding,andk 1186 MOORE AND G•YSON: CATCHMENT PARTITIONING AND RUNOFF PREDICTION TABLE 1. Soil/Vegetation HydraulicPropertiesand AdoptedModel Parametersfor the R-5 Catchment at Chickasha, Oklahoma Manning's n SoilType K,mm h-1 Hc, mm Os,m3 m-3 0r, m3 m-3 Renfrow silt loam Grant silt loam 0.10 0.10 15.9 16.0 100 100 0.49 0.49 0.08 0.08 9.0 9.0 Kingfishersilt loam 0.13 13.7 100 0.49 0.08 9.0 From Moore and Foster [1990] and Sharma et al. [1980]. is a constant. Equation (7) is Horton's equation, which Morel-Seytoux [1988b] demonstrated is equivalent to the Green and Ampt equation [Green and Ampt, 1911] if the parameter k is defined in terms of the Green and Ampt parametersK, Os, Oi, and Hc by = - 0 (8) In practice, k is generally a fitted parameter. Equation (7) is usedin the model becauseit gives an explicit solutionfor the average infiltration rate during a time interval, whereas the Green and Ampt equation must be solved implicitly. Equations (6) and (7) assumethat the depth of pondingon the soil surface is negligible. $chmid [1989, p. 7] showed that this assumptionproducesa maximum error of ! 1% (but generally lessthan 5%) in the computedrainfall excessrates and that "the usual neglectof water depth by overland flow modellers appearsjustified." When the rainfall intensity falls below the infiltration capacityat time ta, drainageof water within the soil profile begins. During drainage, the soil water profile is assumedto be rectangular and the relative hydraulic conductivity is related to soil water content by the Brooks and Corey [1964] equation: after the rainfall ir•tensitydropsbelow •, until the computed depth of overland flow at the outlet of the element reducesto zero. This infiltration algorithm is compatationally efficient and has five physicallybased, measureableparametersK, Os, Or, Hc, and/•, and one initial condition Oi. Hortonianoverlandflow is routedbetweenelementsusing essentially the same equations used in the subsurfaceflow and saturation overland flow model described above, i.e., using (2b), (3b) and (4), but with D = 0. Surface runoff continuesafter rainfall ceasesor when the rainfall intensity is less than the prevailing infiltration capacity, and this surfacerunoff can infiltrate in a downslope element if the soil profile of the downslopeelement is not saturated, i.e., where Hortonian overland flow is occurring. This runoff-runon problem is likely to be important where broad sheet flow occurs but relatively unimportant when the flow concentrates in defined rills, channels and streams. Therefore infiltration is assumed to continue and this runoff-runon phenomenonis represented in the model on elements where broad sheet overland flow occurs. The 9.6-ha R-5 catchment at Chickasha, Oklahoma (35øN 98øE), was selectedfor demonstratingthe application of the Hortonian overland flow model. This application has been described in more detail by Moore and Foster [1990]. A topographic map of the catchment showing the partitioning by TAPES-C, which created 299 elements, is presented in Figure 4. About 65% of the area has slopesof less than 3% where •0 is the hydraulicconductivityat a soil water content and the maximum slope is about 8%. The soils on the of 0, Oris the residualsoil water content,and/• = 3 + 2/qo, catchment consist of Renfrow silt loam (51%), Grant silt where qois the pore size index. For any time t > td when loam (43%), and Kingfisher silt loam (6%). The catchmentis i = 0, the hydraulic conductivity can be predicted by in native grass and is under continuous beef cattle grazing. The model parameters were estimated from soil textural properties and measured data reported by Sharma et al. [1980] and are presented in Table 1 [Moore and Foster, 1990]. In averagingthese data within each of the three soil or when0 < i< •by types, the extreme values reported by Sharma et al. [1980] KO "-KOs-Or/ (9) •o F•+I•oe(t-re) (10) were excluded. iC •0 C- 1 (11) where 0•- i e (12) of this our values of saturated tent are higher than those used by Clapp et al. [1983] and Loague and Freeze [1985] in their modeling studies on the R-5 catchment. C= Because hydraulic conductivity are lower and saturated water con- These studies overestimated soil water con- tent and underestimatedpeak runoff rates and runoff volumes, respectively,which is consistentwith an overestimation of saturated hydraulic conductivity. andFa is the infiltrationvolumeat time td, /• is the mean The April 12, 1967,stormevent [Burfordand Clark, 1973] infiltrationvolumebetweentimestd and t, and K0ais the was selectedto demonstratethe capabilitiesof the model. hydraulicconductivityat time te [Morel-Seytoux,1988a]. This value of •0 can be substitutedinto (9) and solvedfor 0 to give the initial soil water content, 0i, for the next infiltration event. The precipitation data were inadequate to allow spatial variationsin rainfall intensitiesto be modeledand so precipitation was assumedto be spatially invariant. About 66 mm of rainfallfell on the catchmentduringthe 5 daysprior to the The rainfallexcessfor eachelementis i e = i - f provided April 12, 1967,event.We crudelyapproximatedthe spatial i > f. infiltrationcontinuesat the infiltrationcapacityrate variation in antecedentsoil water content by assumingthat MOOREANDGRAYSON: CATCHMENT PARTITIONING ANDRUNOFFPREDICTION 0200 0300 25 0400 0500 0600 0700 :. 20 0800 25 i - I ................ Predictedto Pond 1187 -- E Observed Predicted 50 n.- from Pond t E E 0 0200 0:500 0400 Time 0500 0600 0700 0800 (h) Fig. 11. Observedhyetograph,observedrunoff hydrograph,and predictedrunoff hydrographsto the pond and from the pond at the outlet of the R-5 catchmentfor the April 12, 1967, storm event. the initial soil water content of the Renfrow and Grant silt not accurately modeling infiltration into layered soils, (2) loamsoilswas0.31m3 m-3 andthat of the Kingfisher silt impreciseestimatesof the infiltrationparameters,(3) spatial loamsoil,locatedin theregionofthewaterway,was0.35m3 variability of the infiltration parameters within soil types, or m-3. Channelized flow wasassumed to occurthroughany (4) the effectsof the pond at the catchmentoutlet. In the past element on the catchment whenever the total upslope contributing area of the element exceeded4.2 ha. Runoff at the catchmentoutlet is measured by a 3' 1 broad crested V notch weir. Ponding of water behind the weir tendsto attenuatethe there have been difficulties in modeling runoff from this catchment,particularunder low flow conditions.This has, in part, been attributed to problems created by the pond at the outlet of the catchment (D. Woolhiser, personal communi- runoff and affects the time to commencement of runoff at the cation, 1990). One of the strengthsof the modeling approach proposed here is its ability to predict the overland flow hydrograph and the flow velocity and depth at any point in complex threedimensionalterrain. To demonstratethis capability, Figure 12 presentsthe predicted hydrographsat three points, A, B, and C, shown in Figure 4 and at the inlet of the pond. nication, 1990). The predicted runoff hydrographsat the outlet of the Location A is in a convergingsection of the catchment with watershed(frompond)and at the inlet of the pond(to pond) channelized flow, B is in a converging zone with shallow are shownin Figure 11. Predictedrunoff initiationbegins sheet flow, and C is in a diverging zone with shallow sheet later,predictedtimeto peakoccursat the sametime,andthe flow. Location C has a small specific catchment area compeakdischargeis about9.2% higherthanthe observed.The pared to either A or B, and the predicted runoff'hydrograph hydrographrecessionis modeledvery well, suggesting that at C respondsquickly to rainfall excess, rising rapidly and the overland flow routing methodologydevelopedin the receding rapidly once rainfall excess ceases. Downslope, modelis adequate.The initial ratesof rise (upto discharges moving progressivelyfrom C to B to A and finally to the of about6 mm h-]) of the observedandpredictedhydro- catchment outlet, the predicted hydrograph peak is delayed graphsare similar, but the observedrise occursearlier. in time and the peak runoff rate in millimeters per hour Figure11 showsthat predictedrunoffbeginsenteringthe decreases. Figure 12 also shows that in a divergent flow pond at 0239, but measureablepredictedrunoff (0.06 mm region, such as C, the predicted hydrograph recession is h-1) fromtheponddoesnotbeginuntilabout0306.Mea- steeperthan in a region of flow convergence, such as B. A surable observed runoff(0.06mmh-1) doesnotbeginuntil shadedisometricdiagram of the predicted flow velocities at 0255.Overall,the predictedhydrograph riseis steeperthan the time of the peak runoff is also presented in Figure 13. In this application of the model we only attempted to the observed,and this couldbe causedby (1) simpleinfiltrationequationssuchas the GreenandAmptandHorton represent large-scale spatial variability in soil properties; equations that assumeverticallyhomogeneous soilprofiles that is, the variation between the three soil types. We did not catchment outlet. The predicted runoff to the pond was routed through the pond using a numericalsolutionof the levelpool storageroutingmethod[Chow et al., 1988]on the basisof stage-storageand stage-discharge relationshipsprovided by D. Goodrich and D. Woolhiser (personalcommu- 1188 MOORE AND GRAYSON: CATCHMENT PARTITIONING AND RUNOFF PREDICTION i '• .............. Catchment Outlet 25 '1 ,:,!! ', •5 i' 5 0 0200 0•00 0400 0500 J -- 0600 0700 Time (h) Fig. 12. Predictedrunoffhydrographsat the catchmentoutlet (to the pond)and at locationsA, B and C (see Figure 4) for the April 12, 1967, storm event. attempt to represent the "within soil type" variability, althoughthe model is capable of representingthis detail of catchment. The two models based on TAPES-C provide information about the responseof any point on the catchvariability, were it known. Overall, the shapesof the pre- ment and so can be tested to ensure that the "right" results dictedand observedhydrographsare similar and in reason- are being obtained for the "right" reasons. able agreement. To demonstratethe effect that significant differencesin soil propertiescould have on the distribution CONCLUSIONS of overlandflow in the landscape,we increasedthe resaturated hydraulicconductivity,•c,of the Renfrow silt loam by Terrain analysisis becomingincreasinglyimportant in the 25%to 20 mm h-• andleft all othermodelparametershydrological,geomorphologicaland ecologicalsciencesfor unchanged. Increasing •c for the Renfrow silt loam soil examining the spatial relationships between processesocdecreased both the peak runoff rate and the total runoff curring on the land surface and in the shallow, subsurface volumesignificantlyat the catchmentoutlet. The soil type regime. A variety of terrain analysis methods are available upslopefrom C is almostentirely Renfrow silt loam, and the includinggrid-based,triangular irregular network-basedand initiationof runoffat C wasdelayedby about17min, andthe contour or vector-based methods, and each has advantages peak runoff rate was reduced 15%. and disadvantagesin different circumstances. The contour In this analysis the effect of changingthe infiltration or vector-basedterrain analysis method, TAPES-C, that is parametersof one soiltype was shownto significantlyaffect presented here permits a realistic representation of the the simulated results. It should be recognizedthat the three-dimensionalnature of natural landscapesfor dynamic variability of real soilswithin a soil type will also have a hydrologic modeling under the constraints of maintaining majorinfluenceon the catchmentresponse.TAPES-C hasa physical realism and reducing the computationalrequirestructurethat definessmallelementsin whichit maybe valid ments. The major advantage of TAPES-C over alternative to assume constant parameter values. However, this cannot be fully exploited due to either a lack of data or absenceof methodsis that it allows the governingoverlandflow and shallowsubsurfaceflow equationsto be reducedto a series a suitablemethodof spatialaveraging.Until theseproblems of coupledone-dimensionalequations,thus permittingcomare overcome, care must be taken in interpreting results fromthesemodelsasthe parametervaluesmayappearto be a physicalpropertyof a particularsoilat a particularpoint but in fact representa regionalresponse[Beven,1989].Some knowledge of thecatchment response asa whole,notjust at the catchmentoutlet,is desirableto ensurethat the processesbeing simulatedare in fact those occurringin the putationally efficient solutions to be obtained. Results from the two simple models demonstrate that terrain-basedmodelsof hydrologicprocessescan be power- ful tools for use in computationalhydrology. The results show that the subdivision of a catchment into elements using the "stream tube" approachis a viable methodof considering the effectsof spatialvariability in dynamichydrologic MOORE ANDGRAYSON: CATCHMENT PARTITIONING ANDRUNOFF PREDICTION / Flow Velocity(rns-•) ABOVE 0.300 .... :•: '? :? 0.100 - 0.300 !,•:l;l:•:i:;,,0.070 - 0.100 i•'•:.:•'•':,.i•0.050 - 0.07'0 ....... ii!i!ii!ii!iiii BELOW 0.050 ! 189 Acl,nowledgments. Publishedas Paper 17,799 of the Scientific JournalSeriesof the MinnesotaAgriculturalExperimentStationon researchconductedunder Minnesota Agricultural Experiment Station ProjectsMIN-12-055, MIN-12-30, and MIN-120056 (Regional ProjectS-211).This studywas supportedin part by the Minnesota Agricultural Experiment Station, grants SES-8912042 and SES8912938 from the National Science Foundation and the Water ResearchFoundationof Australia. Travel funds provided to the authorsby theCentrefor Environmental AppliedHydrologyandthe MinnesotaAgriculturalExperimentStationare also gratefullyacknowledged.The authors wish to thank David Goodrich and David Woolhiser of the USDA-ARS, Tucson, Arizona, for the stagestorageand stage-dischargecharacteristicsof the pond and Arlin Nicks of the USDA-ARS, Durant, Oklahoma, for the soil and vegetationdata for the R-5 catchmentat Chickasha,Oklahoma. The commentsof three anonymousreviewers were appreciated. REFERENCES Abbott, M. B., J. C. Bathurst,J. A. Cunge,P. E. O'Connell, and J. Rasmussen,An introductionto the European System-System HydrologiqueEuropeen, "SHE," 2, Structureof a physicallybased,distributedmodelingsystem,J. Hydrol., 87, 61-77, 1986. Anderson,M. G., andT. P. Burt, A laboratorymodelto investigate the soil moistureconditionson a drainageslope,J. Hydroi., 33, 383-390, 1977. Band, L. E., Topographicpartition of watershedswith digital elevation models, Water Resour. Res., 22, 15-24, 1986. Beasley, D. B., L. F. Huggins, and E. J. Monke, ANSWERS: A model for watershedplanning, Trans. ASAE, 24, 938-944, 1980. Betson,R. P., What is watershedrunoff?.,J. Geophys.Res., 69, 1541-1551, 1964. Fig. 13. Isometric projectionof the predicteddistributedflow velocitiesat the time of peak runofffor the April 12, 1967,storm Beven, K., Kinematic subsurface stormflow, Water Resour. Res., event on catchment Beven, K., On subsurfacestormflow:An analysisof responsetimes, R-5. 17, 1419-1424, 1981. Hydrol. Sci. J., 4, 505-521, 1982a. Beven, K., On subsurfacestormflow:Predictionswith simplekinematic theory for saturated and unsaturated flows, Water Resour. Res., 18, 1627-1633. 1982b. modeling.In particular, these modelsprovide a meansof determiningrunoffrates, flow velocities,flow depths,and Beven, K., Changingideas in hydrology•The case of physicallybased models, J. Hydroi., 105, 157-172, 1989. Beven, K. J., and M. J. Kirkby, A physically-basedvariable contributingarea modelof basinhydrology,Hvdroi. Sci. Bull., 24, saturatedareas throughout complex three-dimensionalland43-69, 1979. scapes.Estimates of these quantitiesare essentialif accurate predictionsof erosion,transport,and depositionof sediment Boughton,W. C., Evaluating partial areas of watershed runoff, J. Irrig. Drain Eng., 113, 356-366, 1987. andcontaminantsin the landscapeare to be achieved. Brakensiek, D. L., Kinematic flood routing, Trans. ASAE, 10, Conceptually,the assumptionthat the potentialgradients drivingthe movementof water in the landscapecan be approximatedby the land surfacegradientseemsreasonable for overlandflow modeling.However, the resultsfrom the subsurface flow model suggestthat this may not be appropriatefor modelingshallow subsurfaceflow in landscapes with high effectivehydraulic conductivities.This question will be pursuedin more depthin a subsequent paper. The contour-basedform of distributedhydrologicmodelingallowsthe hydrologistto examinethe effectsof spatial variabilityof soil and vegetationpropertieson runoffgenerationand contaminanttransportprocesses. In addition,the hydrologicresponseof any point in the catchmentcan be Cayley,A., On contourand slopelines,Philos.Mug., 18, 264-268, computed so observation and measurement of the distrib- Dunne, T., and R. D. Black, Partial-area contributions utedcatchment behavior can ensurethat the outflow hydrographis beingproducedby the actualprocessesoccurringon 340-343, 1967. Brooks, R. H., and A. T. Corey, Hydraulic propertiesof porous media, Hydrol. Pap. 3, 27 pp., Colo. State Univ., Fort Collins, 1964. Burford,J. B., and J. M. Clark, Hydrologicdata for experimental agricultural watersheds in the United States 1967, Misc. Publ. 1262,514 pp., USDA-Agric. Res. Serv., Washington,D.C., 1973. 1859. Chow, V. T., D. R. Maidment, and L. W. Mays, Applied Hydroiogy, pp. 245-252, McGraw-Hill, New York, 1988. Clapp,R. B., G. M. Hornberger,andB. J. Cosby,Estimatingspatial variability in soil moisture with a simple dynamic model, Water Resour. Res., 19, 739-745, 1983. Dunne,T., Relationof fieldstudiesand modelingin the predictionof storm runoff, J. Hydrol., 65, 25-48, 1983. to storm runoffin a smallNew Englandwatershed,Water Resour.Res., 6, 1296-1311, 1970. P. W., The hydrologic function of forest soils in watershed thecatchment,i.e., that the modelgivesthe rightresultsfor Fletcher, management, J. For., 50, 359-362, 1952. therightreasons.As discussed by Klemes[1986]andBeven Freeze, R. A., Role of subsurfaceflow in generatingsurface runoff, [1989],thisis vitalfor the development of anytrulyprocess- 2, Upstream source areas, Water Resour. Res., 8, 1272-1283, 1972. basedhydrologicmodel.Althoughnot demonstrated here, R. A., Stormflowgeneration,Rev. Geophys., 12, 627--647, the effectsof spatiallyand temporallyvaryingrainfallpat- Freeze, 1974. terns on catchment hydrology can also be examined rela- Grayson, R. B., Terrain based hydrologicmodelingfor erosion tively easily usingthe detailedcatchmentdiscretization studies, Ph.D. thesis, Dep. Civ. and Agric. Eng., Univ. of methodprovidedby TAPES-C. Melbourne, Parkville, Victoria, Australia, 1990. 1190 MOORE ANDGRAYSON: CATCHMENT PARTITIONING ANDRUNOFF PREDICTION Green,N.H., andC. A. Ampt,Flowof airandwaterthrough soils, J. Agric. Sci., 4, 1-24, 1911. Hewlett,J. D., andA. R. Hibbert,Factors affecting theresponse of by H. J. Morel-Seytouxand D. G. DeCoursey,pp. 248-259, Colorado State University, Fort Collins, 1988b. Mosley,M.P., Streamflowgenerationin a forestedwatershed,New small watershedsto precipitationin humid areas, in Forest Zealand, Water Resour. Res., 15, 795-806, 1979. Pergamon,New York, 1967. relationsto soil and topographicproperties,J. Hydrol., 53, Hydrology, editedbyW. E. Sopper andH. W. Lull,pp.275-290, O'Loughlin, E. M., Saturationregionsin catchmentsand their 229-246, 1981. Hewlett,J. D., andW. L. Nutter,The varyingsourceareasof streamflow fromuplandbasins,in Symposium on Interdiscipli- O'Loughlin,E. M., Predictionof surfacesaturationzonesin natural naryAspectsof Watershed Management, pp. 65-83,American catchmentsby topographicanalysis, Water Resour. Res., 22, 229-246, 1986. Societyof Civil Engineers, Bozeman,Mont., 1970. C.A.,Watershed flood routing using •istributed parameHorton,R. E., The roleof infiltration in the hydrologic cycle,Eos Onstad, Trans. AGU, 14, 446-460, 1933. Jones,J. A. A., Thenatureof soilpiping:A reviewof research, in BritishGeomorphological Research GroupResearch Monograph No. 3,301 pp., Geo Books,Norwich, 1981. Kibler, D. F., and D. A. Woolhiser,The kinematiccascadeas a hydrologic model,Hydrol.Pap. 39, 27 pp., Colo.StateUniv., ters, in Floods and Droughts, edited by E. S. Schulz, V. A. Koelzer, and K. Mahmood, pp. 418-428, Water ResourcesPublications, Fort Collins, Colo., 1973. Onstad,C. A., and D. L. Brakensiek,Watershed simulationby the streampath analogy,Water Resour.Res., 4, 965-971, 1968. Quimpo,R. G., Spatialheterogeneity andmodelsof surfacerunoff, Fort Collins, 1970. J. Hydrol., 68, 19-28, 1984. Water Resour. Res., 22, 177S-188S, 1986. treatedas independentof flow depth?,J. Hydrol., 107, 1-8, 1989. Klemes,V., Delettantismin hydrology:Transitionor destiny?, Schmid,B. H., On overland flow modeling: Can rainfall excessbe Knapp,B. J., Infiltrationand storageof soilwater,in Hillslope Sharma,M. L., G. A. Gander, and C. G. Hunt, Spatialvariabilityof infiltration in a watershed, J. Hydrol., 45, 101-122, 1980. Hydrology,editedby M. J. Kirkby,pp.43-72,JohnWiley, New Sivapalan,K., K. Beven,andE. F. Wood, On hydrologicsimilarity, York, 1978. Knisel,W. G. (Ed.), CREAMS: A fieldscalemodelfor chemicals, 2, A scalemodelof stormrunoff production, Water Resour.Res., 23, 2266-2278, 1987. runoffand erosionfrom agriculturalmanagement systems,Conserv. Res. Rep. 26, 643 pp., USDA-Sci.and Educ. Admin., Skaggs,R. W., andR. Khaleel,Infiltration,in HydrologicModeling of Small Watersheds,ASAE Monogr., vol. 5, editedby C. T. Washington, D.C., 1980. Haan, H. P. Johnson, and D. L. Brakensiek, pp. 121-166, Kozak, M., Determinationof the runoffhydrographon a determinAmerican Society of Agricultural Engineers, St. Joseph,Mich., istic basisusinga digitalcomputer,IAHS Publ., 80, 138-151, 1982. 1968. Lane, L. J., E. D. Shirley, and V. P. Singh,Modelingerosionon Sloan,P. G., and I. D. Moore, Modelling subsurfacestormflowon steeply slopingforested watersheds, Water Resour. Res., 20, hillslopes,in ModellingGeomorphological Systems,editedby 1815-1822, 1984. M. G. Anderson,pp. 287-308,JohnWiley, New York, 1988. Lee, M. T., andJ. W. Dellcur, A variablesourceareamodelof the Smith, R. E., and R. H. B. Hebbert, Mathematicalsimulationof interdependentsurfaceand subsurfacehydrologicprocesses,Warainfallrunoffprocessbasedon the watershedstreamnetwork, ter Resour. Res., 19, 987-1001, 1983. Water Resour. Res., 12, 1029-1036, 1976. Loague,K. M., Streamflowgeneration:Equallylikely realizations Smith, R. E., and D. A. Woolhiser,Overland flow on an infiltrating surface, Water Resour. Res., 7, 899-913, 1971. of a singlestochastic process,in ModelingAgricultural,Forest, and RangelandHydrology,pp. 420-426,AmericanSocietyof Speight,J. G., A parametricapproachto landformregions,Spec. Publ. Inst. Br. Geogr., 7, 213-230, 1974. AgriculturalEngineers,St. Joseph,Mich., 1988. Loague,K. M., andR. A. Freeze,A comparison of rainfall-runoff Speight,J. G., The role of topographyin controllingthroughflow generation:A discussion, Earth Surf. Processes Landforms,5, modelingtechniqueson smalluplandcatchmerits,WaterResour. Res., 21,229-248, 187-19 I, 1980. 1985. Mark, D. M., Conceptsof "data structure" for digital terrain Sunada,K., and T. F. Hong, Effectsof slopeconditionson direct runoffcharacteristics by the interflowandoverlandflow model,J. models,paperpresentedat DigitalTerrainModels(DTM) SymHydrol., 102, 323-334, 1988. posium,Am. Soc.of Photogramm. St. Louis,Mo., May 1978. Mein, R. G., and C. L. Larson,Modelinginfiltrationduringa steady Takasao,T., and M. Shiiba, Incorporationof the effect of concentrationof flow into the kinematicwave equationsand its applicarain, Water Resour. Res., 9, 384-394, 1973. tionto runoffsystemlumping,J. Hydrol., 102,301-322,1988. Moore, I.D., and G. J. Burch, Sedimenttransportcapacityof sheet and rill flow: Application of unit stream power theory, Water Takasao,T., M. Shiiba,and H. Kitamura, A distributedrunoff modelandlumping scaleof a fiverbasin,Proc.31stAnnu.Conf. Resour. Res., 22, 1350-1360, 1986. Moore, I. D., and G. R. Foster, Hydraulicsand overlandflow, in ProcessStudiesin Hillslope Hydrology, edited by M. G. Anderson and T. P. Burr, pp. 215-254, JohnWiley, New York, 1990. Moore, I. D., and R. B. Grayson,Hydrologicand digitalterrain modellingusingvectorelevationdata(abstract),Eos Trans.AGU, 70, 1091, 1989. Moore, I.D., S. M. Mackay, P. J. Wallbrink, G. J. Burch, and E. M. O'Loughlin, Hydrologiccharacteristics and modellingof a small forested catchment in south-eastern New South Wales: Pre- loggingcondition,J. Hydrol., 83, 307-335, 1986. Jpn. Soc. Civ. Eng., 2-86, 1-163, 1976. Tanaka,T., The role of subsurface water exfiltrationin soilerosion processes,IAHS Publ., 137, 73-80, 1982. Tanaka,T., M. Yasuhara,H. Sakai,andA. Marui, The Hachioji experimental basinstudy---Storm runoffprocesses andthemechanismof its generation,J. Hydrol., 102, 139-164,1988. Tisdale, T. S., J. M. Hamrick, and S. L. Yu, Kinematicwave analysisof overlandflow usingtopographyfitted coordinates (abstract),Eos Trans. AGU, 67, 271, 1986. Van de Griend,A. A., and E. T. Engman,Partialarea hydrology and remote sensing,J. Hydrol., 81,211-251, 1985. topographicmodelfor hydrologicaland ecologicalapplications, Vieux, B. E., V. F. Bralts, and L. J. Segerlind,Finite element analysis of hydrologic response areasusinggeographic informaEarth Surf. ProcessesLandforms, 13, 305-320, 1988. tionsystems, in ModelingAgricultural, Forest,andRangeland Moore, i.D., R. B. Grayson, and A. R. Ladson, Digital terrain Hydrology, pp.437-446,American Societyof Agricultural Engimodelling:A review of hydrological,geomorphological and bioneers, St. Joseph,Mich., 1988. logicalapplications,Hydrol. IProc.,5, 3-30, 1991. andcontourmapping, J. Hydrol.,25, Morel-Seytoux, H. J., Recipe for simple but physically based Warntz,W., Streamordering 209-227, 1975. modelingof the infiltrationand local runoff process,in ProceedD. R., Measurements of the downslope movement of ingsof the EighthAnnualAGU Front RangeBranchHydrology Weyman, water in a soil, J. Hydrol., 20, 267-288, !973. Days, edited by H. J. Morel-Seytouxand D. G. DeCoursey,pp. Moore, I. D., E. M. O'Loughlin, and G. J. Burch, A contour-based Whipkey,R. Z., Stormrunofffromforestedcatchments by subsur226--247,Colorado State University, Fort Collins, 1988a. face routes, IAHS Publ., 85, 773-779, 1969. Morel-Seytoux, H. J., Equivalencebetweeninfiltrationparameters in Horton and Morel-Seytouxformulae,in Proceedings of the EighthAnnual AGU Front RangeBranchHydrologyDays, edited Woods,R. A., andR. P. Ibbitt,Analyticalsolution for kinematic flowover an infiltratingplane,in ModelingAgricultural,Forest, MOORE ANDGRAYSON: CATCHMENT PARTITIONING ANDRUNOFF PREDICTION andRangeland Hydrology, pp.61-71,American Society of Ag- riculturalEngineers,St. Joseph,Mich., 1988. I. D. Moore, Centre for Resource and Environmental Studies, AustralianNational University, G. P.O. Box 4, Canberra, ACT Young, R. A., C. A. Ohstad, D. D. Bosch,andW. P. Anderson, 2601, Australia. AGNPS:A nonpoint-source pollution modelforevaluating agricultuml watersheds, J. SoilWaterConserr.,44, 168--173, 1989. R.B.Grayson, Department ofCivilandAgricultural Engineering, Centre for Environmental AppliedHydrology, University of Mel- bourne,Parkville, Victoria 3052, Australia. 1191 (Received January 18, 1990; revised September 28, 1990; accepted January 8, 1991.)