key: cord-0009674-plxq9e7x authors: Moore, I. D.; O'Loughlin, E. M.; Burch, G. J. title: A contour‐based topographic model for hydrological and ecological applications date: 2006-07-18 journal: Earth Surf Process Landf DOI: 10.1002/esp.3290130404 sha: 69f4b99eeb03310c3761b6ff7b1d1119680e7022 doc_id: 9674 cord_uid: plxq9e7x A digital model for discretizing three‐dimensional terrain into small irregularly shaped polygons or elements based on contour lines and their orthogonals is described. From this subdivision the model estimates a number of topographic attributes for each element including the total upslope contributing area, element area, slope, and aspect. This form of discretization of a catchment produces natural units for problems involving water flow as either a surface or subsurface flow phenomenon. The model therefore has wide potential application for representing the three‐dimensionality of natural terrain and water flow processes in the fields of hydrology, sedimentology, and geomorphology. Three example applications are presented and discussed. They are the prediction of zones of surface saturation, the prediction of the distribution of potential daily solar radiation, and the prediction of zones of erosion and deposition in a catchment. Computer based modelling of the hydrology, sedimentology, and geomorphology of real three-dimensional catchments requires overland flow and/or subsurface flow, which are the prime mechanisms for solute and sediment transport, to be modelled throughout the catchment. The topographic variables required to model these pbenomena include: (1) upslope contributing area; (2) local slope of the potentiometric surface (which can often be approximated by the slope of the land surface); and ( 3 ) local aspect (which together with the local slope determines the potential solar radiation received by any point on a catchment, and hence potential evapotranspiration and snowmelt). As these variables can vary greatly over a catchment, it is important that estimates of them be obtained at locations and at scales that can reflect both the local and integrated effects of topography on runoff generation and catchment soil water status. Digital terrain or digital elevation mod&k (DTMs or DEMs) are the most common methods used for automatically extracting these topographic variables from raster elevation data. Manual methods of extraction, such as that used by Beven and Kirkby (1979) , are time consuming and intractable on all but relativelysimplecatchments. Band (1986) and O'Callaghan and Mark (1984) describe theapplication of DEMs for determining drainage divides and stream networks. However, their discretization of catchments is too coarse for detailed modelling of runoff processes in many applications. The most widely used data structures for DTMs and DEMs consist of grid networks. For example, the methods of determining shape and aspect developed by Sharpnack and Akin (1969) and Travis el al. (1975) 01 97-9337/88/040305 -1 6$08.00 0 1988 by John Wiley & Sons, Ltd and Clerici's (19801 and Mulla's (1986) method of deriving slope maps are based on grid cells. Heerdegen and Beran (1982) fitted five-parameter polynomials to grid elevation data derived from contour maps to compute both horizontal and vertical curvature characteristics of catchments. Armstrong (1976) , Ahnert (1976) and Hirano (1976) used grid networks for determining and inputting the topographic variables into their simple three-dimensional slope models of landform development that included submodels of the overland flow process. The ANSWERS model (Beasley et a/., 1980) is one of the few hydrology/erosion/water quality models capable of representing the three-dimensionality of landscapes, but it too is based on a grid network. Most grid cell methods developed to date for hydrologic application have been too simplistic, even when used at very fine scales. Their principal drawbacks are that: (1) they do not allow water flow from one cell to be split, leading to significant error in divergent areas; and (2) the directions of flow trajectories are matched only crudely by transitions from one grid cell to another, even in very simple planar regions. Mark (1978) noted that grid structures for spatially partitioning topographic data are not appropriate for many geomorphological applications, and in particular for digital terrain modelling. He stated that 'the chief source of this structure should be the phenomena in question, and not problems, data, or machine considerations, as is often the case'. Natural units into which a catchment should be subdivided to represent the phenomenon of waterflow are polygons formed by equipotential lines and their orthogonals, streamlines. For many water-flow processes occurring on three-dimensional catchments these can be approximated by contour lines and their orthogonals (flow trajectories) that define the boundaries of upslope drainage areas. The disadvantage of this approach is that one needs at least an order of magnitude more points in contour line form than in regular grid form to adequately describe an elevation surface. Also, it is computationally slower than the grid cell approach. This paper describes a contour-based model that partitions three-dimensional catchments into natural units consisting of irregularly shaped polygons and estimates upslope contributing area, slope, and aspect for each unit. The distributed topographic variables calculated by the model have been used in several applications which include the identification of zones of saturation (O'Loughlin, 1986 ) and zones of erosion and deposition (Moore and Burch, 1986b) and the estimation of potential daily solar radiation in real three-dimensional catchments. Examples of these applications are presented. An idealized topographic map of a small catchment is presented in Figure 1 . The model was designed to calculate, for each section of a contour (an example of which is the section defined by the line joining j , j + 1 in Figure l), the upslope trajectories or streamlines from the segment end-points to high point(s) on the catchment boundary ( j , 0 and j + 1, 0 in Figure 1 ). For each section of contour the following topographic attributes are estimated: (1) the total upslope contributing area, A j , bounded by the contour section j , j + 1 and the pair of adjacent trajectories, j , 0 and j + 1, 0; (2) the area of the contour element bounded by adjacent contours and adjacent trajectories, A , (e.g. polygon defined by the points i, i+ 1, j + 1, j in Figure 1 ); (3) the widths of the contour sections defining A , (bj and &in Figure 1) ; (4) the average local slope orthogonal to the contour along each contour segment (i.e. at Pj and Q j in Figure 1 ); and (5) the x, y , z coordinates of P, and R , (which is the centroid of the element defining A,) in Figure 1 . The model approximates contours by short straight line segments and trajectories between adjacent contours as straight lines. These approximations were made to simplify the model so that no iterative techniques were used in determining the trajectories, thus reducing the computer processing time, and to minimize storage requirements for arrays. These approximations cause little error in the estimation of the trajectories when contours are close together, but can produce substantial error when they are not. This model has evolved from a topographic model originally developed by O'Loughlin (in press). The model calculations are performed in two programs; PREPROC and TOPO. PREPROC preprocesses the input digitized contour information, while TOPO is the program in which the trajectories and topographic attributes, described above, are calculated. Each of these programs is briefly described below. The model requires sets of x, y coordinates for each contour of elevation z, the x, y, z coordinates of all high points on the catchment boundary, and the x, y coordinates of points defining the catchment boundary. It can accept up to 1500 sets of x, y coordinates for each contour line, up to 80 contour lines, up to 10 high points (on the catchment boundary), and up to 40 points defining the catchment boundary. These limits can be increased by modifying the dimension statements in the program. Any high points on the catchment boundary must also be included within the 40 sets of boundary points. The boundary points are input into the model in program TOPO. The contour and high point data for input into the model can be generated in two ways; either by digitizing an existing topographic map, or by applyinga DEM to generate contours in digital form from a regular grid of raster elevation data. If a topographic map of the catchment is already available then each contour on the map can be sequentially digitized using a flat-bed digitizer, beginning with the contour of lowest elevation. First, two points on a reference baseline of known azimuth are digitized. These are used to transform the raw digitized coordinate pairs (x, y coordinates) into a set of true eastings and northings. Second, the elevation of each contour is keyed-in (e.g. CON40, for aicontour of elevation 40 m) and the contour is then digitized, producing a string of up to 1500 sets of irregularly spaced coordinates for each contour. A high point is input by keying-in the prefix HPT, followed by its elevation and then digitizing the (x, y) coordinates of the point. A high point is entered into the input file immediately after digitizing the contour segment above which it lies. If the elevation data are in the form of randomly distributed x, y, z coordinates then a DEM consisting o f a regular grid of elevations can be fitted to these data. Digital contours and high point data for input to the model can then be calculated from this regular grid. The authors have used programs developed by Hutchinson (1981 Hutchinson ( , 1984 to do this. Program PREPROC is a preprocessing program that transforms the raw digitized x, ycoordinates and high points into a form that allows the topographic analysis to be performed. The raw coordinates and high points are transformed into a set of true eastings and northings and rescaled to lie within the bounds of 0 to 1000 in both the x and y directions. If the raw coordinates were generated by digitizing a topographic map then they would be irregularly spaced along a contour. The program linearly interpolates a new string of regularly spaced x, y coordinates defining each contour. The interpolation interval is determined by the user input parameter RINT. The results of such an interpolation is demonstrated in Figure 2 , which is an enlarged segment of a contour. In this figure the original rescaled coordinates defining the contour are shown as solid circles and the interpolated points are the open circles. Users have the option of smoothing abrupt irregularities in the contours by passing the coordinates through a two-pass moving average filter as well as the option of plotting the contours. An example of this plotted output is presented in Figure 3 . The x, y coordinates of all points defining the contours are stored as two contiguous coordinate vectors (one for the x and the other for the y coordinates), beginning with the first point on the lowest contour and ending with the last point on the highest contour. The location of the first point on each contour within these vectors is also stored in vector form. If one or more gaps (greater than a specified size and where a boundary intersects the contour at two or more points) occur in a contour line then it is divided into a number of contour segments and the location of the first point of each segment and the number of segments for each contour are stored in two vectors. This method of storing the contour information permits rapid and efficient retrieval of data during the execution of program TOPO. The trajectories and topographic attributes of a catchment are computed in program TOPO. The catchment boundaries are input as up to 40 sets of x, y coordinates joined by straight line segments (see Figure 4 ). All contour points lying outside the catchment's boundary are excluded from the calculations. These points are determined using Hall's (1975) point in polygon routine by taking the signed summation of the angles subtended at the given point by each segment of the closed catchment boundary. If the angles sum to f 360 then the point is inside the boundary, and if they sum to 0" then the point is outside the boundary. Using this test the first and last points on each contour segment lying on or within the catchment's boundary are identified. Contours are digitized in a consistent direction so that model calculations proceed in the same direction for all contours. The model takes the first point on each contour segment lying on or within the boundary and identifies this as the first trajectory start point. As the calculations proceed each contour segment is divided into small sections of length b, (see Figure I) , determined by the user input parameter PINT, thus identifying the Figure4. An example ofthecomputer generated output from programTOP0 showing the high points (HP)and thecatchment boundary nodes (solid circles). The computed upslope trajectories originating from only three contour lines are shown tor simplicity coordinates of the starting points for all trajectory calculations (e.g. the pointsjandj + 1 in Figure 1 ). PINT is a multiple of RINT (input in program PREPROC) and good results have been obtained by the authors when PINT/RINT > 10. The x,y coordinates of the mid-point of the section, Pi (see Figure l) , are also calculated. In Figure 1 subscripts j and j + 1 etc., refer to a given contour line for which the topographic attributes are being estimated and i and i + 1 etc., refer to its upslope contour. The points Q j and P j are the points midway betweenj a n d j + 1, and i and i + 1, respectively, whereas R j is the point defining the centroid of the elemental area Ar,. Each trajectory point is also a contour point and its address within the vectors containing ?he x, y coordinates of the contours is determined by the program. The points on the contour were interpolated to a regular spacing (RINT) in program PREPROC, so that bi, the distance along the uphill contour between adjacent trajectories, is RINT times the difference between the addresses of the trajectory points i and i + 1. Linsley et al. (1949) define the aspect of a slope as 'the compass direction normal to the slope contour and pointing downslope'. In the model the aspect of the jth section of the contour (i.e. at P j ) is estimated by the orthogonal to the direction of the straight line connecting the point j to the point j + 1 (approximately the tangent of the contour line at P j ) , in the downslope direction. It is initially computed as an angle (0 to 2n radians) anticlockwise from east, and is then adjusted to give a true azimuth (clockwise from north). Upslope trajectories from a point are computed using two criteria: minimum distance, or orthogonals. The two criteria are used in an attempt to overcome the error caused by using straight line segments between adjacent contours to define the streamlines or trajectories. The error increases as the contour lines get further apart. The application of the two criteria is determined based on the curvature of the contour line, and is controlled by a user-input parameter. The minimum distance criterion uses the minimum distance between adjacent contours to define trajectories and is applied to ridge areas where downslope flow diverges. In valleys and where flow converges trajectories are estimated by computing the orthogonal to the point of interest and projecting it to the upslope contour. In practice trajectories are streamlines that are always orthogonal to equipotential lines. A contour line is a line of equal potential energy for most terrain analyses. However, for most surface and subsurface water flow phenomena this is an approximation because the energy grade line in the downslope direction is not necessarily parallel to the land surface. The technique for determining the orthogonal to thejth trajectory start point (see Figure 1 ) is similar to that described above for computing the aspect at P j . The differences are that the orthogonal is in the upslope direction and the points j k h j / 2 are used to compute the approximate direction of the tangent to the contour at j. Calculations for the jth trajectory (see Figure 1 ) begin at the trajectory start point, j. The direction of the orthogonal to j is calculated and then the string of points on the uphill contour within a user determined maximum distance fromj are examined. The distance and direction (0 to 2n radians) betweenj and each point are computed. Points satisfying the two criteria are then selected by the model. If no point satisfying the orthogonal criterion is found, then the point determined by the minimum distance criterion is used. The curvature test is then applied and the correct uphill trajectory point is selected (i in Figure 1 ). The procedure is repeated for successive trajectory points until a high point or ridge, defined as the highest contour (if no high point exists), is reached. An example of trajectories computed by the model for three different contour lines of the catchment shown in Figure 3 is presented in Figure 4 . The model computes two average slope term:'for each elemental area, S1 and S2 (m m-'), in which and where Ah is the contour interval, S1 is the average slope of two adjacent trajectories between a contour section and its first upslope contour, and S2 is the average slope between P j and Q j . On planar slopes and in ridge areas S1 = S2, but if the elemental area straddles a valley S1 # S2. In such cases S1 is an approximation of the slope of the land draining to a perennial or intermittent waterway, while S2 approximates the slope of the waterway itself. The model also estimates the average slope of the catchment between Q j and the first contour above Q j , adopting the same method used to calculate S2. The final topographic attributes calculated by the model in program TOPO are the two area terms A,j, the area of the contour element bounded by adjacent contours and adjacent trajectories, and A j , the upslope contributing area (see Figure 1 ). The area enclosed by a polygon with n vertices at (xk, yk) is: A,j and A j are calculated by applying equation (3) to the points making up the boundaries of the two areas. If the two adjacent trajectories defining A j terminate at two different high points then the (xk, yk) points defining the area must also include the nodes along the catchment boundary between the high points. It is for this reason that high points must also be catchment boundary nodes. The present version of the model can not handle high points that lie within the catchment's boundary, although work is underway to overcome this limitation. The Geebung Creek catchment is used here to demonstrate the ability of the model to predict distributed topographic attributes of a catchment and is also used in the example applications that follow. The catchment is 79.6 ha in area and is located about 30 km inland from the southeast coast of New South Wales, Australia, in the 51 300 ha Yambulla State Forest (37"18'S and 149"40E). It has an easterly aspect (90) and a relief of about 174 m. Soils are coarse textured and are generally less than one metre deep in upslope areas. They are somewhat deeper in lower-slope areas and many rock outcroppings occur throughout the catchment. The soils are highly to very highly erodible and have a low nutrient status. The vegetation consists of dry sclerophyll forest having a tall open structure. The understorey is generally sparse except for localized moderately dense stands of Casuarina littoralis on ridges and dense thickets of Melaleuca squarrosa, Banksia serrata and associated species and a variety of grasses around drainage lines and other wet zones (Mackay and Cornish, 1982; Moore et at., 1986) . Average annual rainfall is in excess of 900 mm per annum. Details of the rainfall-runoff response of the catchment are given by Moore et al. (1986) . A 1 : 5000 scale, 5 m contour interval, topographic map of the catchment (see Figure 3 ) produced from aerial photographs, was digitized and formed the primary input to the model. During the digitizing process six highpoints were identified, as shown in Figure 3 . Eighteen coordinate pairs defining the catchment's boundary (see Figure 4) Three-dimensional projections of the catchment overlaid with each of the three computed topographic attributes are presented in Figures 5, 6 and 7 . These projections were produced by UNIRAS, a set of FORTRAN programs for generating computer graphics. UNIRAS generates a three-dimensional surface with an overlay of colours (or a grey scale) that can be chosen to represent the variation of a location dependent function (in this case one of the three topogra;hic attributes). The software interpolates the topographic data (x,y,z data) to a regular grid and calculates smoothed values of the function over the same region. Figure 5 shows the contributing area (in hectares) above each 25.4 m long contour segment and increasingly darker shadings reflect increasingly larger contributing areas. As expected, the smallest contributing areas are located on ridges (where flow diverges) and topslope sections of the catchment and the largest areas are in the valleys and drainage lines (where flow.converges). The drairiage lines are readily identifiable in this figure as the regions with darker shadings. The relationship between stream order and contributing area for the Geebung Creek catchment is also readily discernible from Figure 5 . First-order streams (or drainage lines) have contributing areas of 1-5 ha and 2nd-order streams seem to form when the contributing area is about 10 ha. Although not shown in this figure, the analysis indicated that the 3rd-order stream is formed when its contributing area exceeds about 50 ha. The Geebung Creek catchment is quite steep, exceeding 40 per cent in some parts, and this is reflected in the large relief of 174 m. The distribution of the various slope classes, computed as S1 using equation ( 2 ) (with slope in per cent), is shown in Figure 6 . The flattest slopes, shown by the lighter shadings, occur on the ridge tops in the upland sections of the catchment and along the valley floor, upslope from the catchment outlet. The steepest slopes, corresponding to the darker shadings, tend to occur in midslope positions between ridges and drainage lines. The dark shading at the catchment outlet reflects an increase in slope of the main stream channel as it falls to meet the Wallagaraugh River. The final topographic attribute computed by the model is aspect, and its distribution across the Geebung Creek catchment is presented in Figure 7 . In this figure aspect is shown in degrees clockwise from north. The catchment has a general easterly aspect, with major portions having northeast to east (45"-90") and east to southeastern (90"-135") aspects. Only small areas have south to north ( 180"-360') aspects, the unshaded areas in Figure 7 , but these are mostly hidden in the figure because of the view angle used to present the results. Therefore, large sections of the catchment would receive the morning sun and would be shaded in the late afternoon. Many aspects of hydrologic behaviour in small catchments are associated with their saturated source area characteristics. These saturated source areas expand and contract as a direct consequence ofthe wetting-up and drying-out ofthe catchment and generate overland flow during rainfall events (Hewlett and Nutter, 1970) . This overland flow appears as rapid runoff in the storm hydrograph. The same areas are often the most sensitive to mechanical disturbance because of low soil strength, and are liable to become salinized if solutes present in the seepage water become sufficiently concentrated at the surface. The location and size of these zones is therefore of direct interest for interpreting a range of catchment hydrologic processes, especially in cases where land use or vegetation management could alter the behaviour of these zones. Beven and Kirkby (1979) , O'Loughlin (198 1,1986) and Zaslavsky and Sinai (1981) observed that variations in wetness within a catchment are explicable imkrrns of the local topography (slope and degree ofconvergence of the hillslope) and the hydraulic properties of the soil profile (uiz. the often observed decrease in hydraulic conductivity with depth below the soil surface). Local saturation at any point in a catchment will occur whenever the drainage flux from upslope exceeds the capacity of the soil profile to transmit the flux. O'Loughlin (I98 1) expressed this criterion as where b is the length of the contour element (see Figure I) , Qb is the local drainage flux across this element, S is the local hillslope gradient (m m-* ), and Tis the soil transmissivity (the depth integrated saturated hydraulic conductivity). This relationship assumes that lateral drainage takes place over an effectively impermeable layer (relative to the lateral water flux) and that the slope of the land surface approximates the water table gradient. By assuming steady-state drainage conditions, O'Loughlin (1986) recast equation (4) in a form that can account for variable drainage fluxes and transmissivities over a catchment. This relationship can be written in a dimensionless form as where q is the drainage flux per unit area (flux density), dA is the area contributing to drainage upslope of a contour element of width b, T and qo are the average catchment transmissivity and drainage flux density, respectively, L is a characteristic length used to make the expression dimensionless (mean hillslope length), D is a drainage index, and W is a wetness index. More detailed descriptions of the derivation of this equation and the assumptions involved are given in O'Loughlin (1986) . The term on the left hand side of equation (5), the drainage index, consists of two dimensionless ratios (q/qo and TIT) and four topographic attributes ( S , b, dA, and L), three of which (S, b, and dA) can be calculated throughout a catchment by the topographic model proposed herein. Thus, by combining equation (5) with the topographic model it is possible to calculate drainage indices, and hence identify zones of surface saturation, in complex three-dimensional terrain. The simplest case that can be considered is to assume that the drainage flux density and transmissivity are uniform throughout the catchment (i.e. q/qo = 1 and T/T = I), in which case equation (5) reduces to so that the drainage index is a function of topographic attributes only. These assumptions are used here simply to demonstrate the application of the method. The predicted locations and sizes of the zones of surface saturation for the Geebung Creek catchment are presented in Figure 8 as a function of the drainage index, D. Increasingly darker shadings in this figure reflect increasingly wetter zones in the catchment for a given average drainage flux density, g o . A given location is predicted to be saturated at the soil surface when the drainage index is greater than the value given in Figure 8 . Observations on the catchment confirm these general predictions (Moore et al., 1986) . Moore et al. (1986) used the predicted wetness indices, together with measured rainfall-runoff data, to estimate the average transmissivity, T, of the Geebung Creek catchment and found good agreement with the values computed from measured hydraulic conductivities obtained from soil cores and in-situ measurements using a well permeameter (on a grid network across the catchment). Furthermore, they used the calculated wetness index versus per cent saturated source area relationship as the basis for a full dynamic simulation of the rainfall-runoff response of the catchment by assuming successive steady-state conditions. The predicted runoff hydrographs showed excellent agreement with the observed. Detailed descriptions of the techniques used are presented by Moore et al. (1986) . These examples illustrate only two ways in which the saturation zone analysis can be used. Further examples include the use of the results by forest managers to optimize the exclusion zones where logging should not be permitted at all and the identification of zones where logging may or may not be permitted, depending on the weather. O'Loughlin (1986) has used the analysis to determine the impacts of selective clearing on the local water balance, and Burch et al. (1987) have used it to differentiate between the effects of topography and soil hydraulic properties on runoff generation in cleared and forested catchments. The effect of topography on the amount of radiation received at the land surface is particularly important in understanding the ecology and the hydrology of catchments. For example, Tajchman and Lacey (1986) , using Budyko's (1974) radiation index of dryness (net radiation/latent heat equivalent of precipitation), showed that biomass production on two small catchments was related to aspect, radiation and wetness. Austin et al. (1983) have also demonstrated that the distribution of some eucalypt species in southeastern Australia is related to altitude, rainfall, and annual radiation index. The potential daily solar radiation, R , is a function of the topographic attributes of slope ( E in degrees) and aspect or azimuth (A'-measured clockwise from north), and the solar declination (6), terrestrial latitude (4), the ratio of the earth-sun distance to its mean (r), the transmission and scattering properties of the atmosphere, and the albedo of the surface (Robinson, 1966; Lee, 1978) . A first and simple assumption for demonstration purposes is to neglect atmospheric and surface effects so that diffuse radiation is ignored. The solar declination and the ratio of the earth-sun distance to its mean are essentially constant over a day so that the potential daily solar radiation can be written as 24 l o r2 R = -(cos @ cos 6) (sin wtlwtlcos w t l ) where lo is the solar constant (4.871 MJ m-' h-' or 1.942 Ly min-I), w is the angular velocity of the earth (412 radians h-'), ti is the sunrise-sunset time from solar noon, as seen by the inclined surface. [coswt; = -tan 4' tan 61, sin 4' = sin ECOS A'cos 4' + cos E sin 4, and the other symbols are as previously defined. The solar declination ranges from 23.5" at the winter solstice to -23.5" at the summer solstice in the southern hemisphere and the ratio of the earth-sun distance to its mean varies in a narrow range from 0.9833 to 1.0167. The topographic variables required to model the distribution of potential solar radiation, R , across a catchment are local slope and aspect, which are properties of each contour segment. Upslope contributing area is not used in this particular application. The distribution of R across the Geebung Creek catchment was calculated at both the winter and summer sdstices (6 = 23.5" and -23.5", respectively), representing the minimum and maximum potential daily solar radiation received by the catchment during a year. The results for the winter solstice are presented in Figure 9 . The potential daily solar radiation on a horizontal plane ( R h ) in the catchment at the winter solstice is 14.1 MJm-'d-' (336 Ly min-'). The R/Rh ratio ranges from a low of 025 on sites with southwesterly aspects to a high of about 1.2 on sloping sites with northern and eastern aspects. At the summer solstice Rh is 44.2 MJ m-' d -' (1058 Ly min-'), more than three-fold the winter solstice value, but the R / R h ratios exhibit less variation, ranging from 0.7 to 1.0. Austin (personal communication: unpublished manuscript) and Binns (1984) have measured the species composition of the overstorey for the Geebung Creek catchment on 0.1 ha circular plots arranged on a grid pattern. A simplified form of these results is presented in Figure 10 , in which the eucalypt species have been divided into four broad groups: (1) legend of Figure 10 also has a high probability of being found at the sites indicated on the figure. Species groups 2, 3, and 4 occur predominantly in the wetter and shadier sites, as shown on Figures 8 and 9 , respectively, while the group 1 species, the E. sieberi, are found mostly in the drier areas receiving higher amounts of radiation. The ridges on the catchment's boundary contain E. sieberi almost exclusively. These ridges are both the driest parts of the catchment, as well as the most exposed, receiving the highest potential daily radiation. On the other hand, the authors have observed that none of the eucalypt species in groups 1 to 4 are found in the wettest drainage lines shown in Figure 8 (the areas with a drainage index above 16.0). These observed distributions may also be affected by fire, which undoubtedly plays a part in the relative abundance among eucalypt species in this forest (Bridges, 1983) . Both radiation and drainage index variations could be expected to influence the susceptibility of parts of the landscape to fire through their effect on fuel moisture content. These qualitative correlations suggest that linking the distributed results from the saturation-zone EUCALYPT SPECIES "4-0 E muellerana (oblrqua cypellocarpa 0 0 0 0 0 0 0 0 or consideniana) analysis with those from the radiation analysis may prove useful in predicting the distribution and location of dift'erent tree and plant species on a catchment. Estimates of erosion rates, the redistribution of soil within the landscape, and sediment loads in runoff are being increasingly required by land and water resource managers to assess the possible impacts of land use and land management practices on water quality, and land productivity and degradation. Recent emphasis in models for predicting sediment transport have concentrated on detailed representation of the physical processes of soil detachment, entrainment, and transport by rainfall and sheet and rill flow at the plot or hillslope scale (for example, Alonso et al., 1981; Rose et a/., 1983; Gilley et al., 1985; and Moore and Burch, 1986a) . These models generally idealize the topography as a two-dimensional hillslope and poorly represent overland flow convergence and divergence. Even though they are quite powerful models this is a major deficiency that makes them incapable of accurately predicting erosion and deposition in naturally complex three-dimensional terrain. However, improved predictions of the sediment fluxes in such complex three-dimensional catchments are possible by combining a soil erosion model with a quantitative three-dimensional description of the hydrology of a catchment, based on a DTM or topographic model. By dividing a catchment into small elements or polygons, the erosion and deposition rates per unit area can be calculated for each element by performing a mass balance of the sediment entering and exiting each element per unit time. The topographic model described in this paper provides a mechanism for dividing catchments into elements compatible with the physics of overland flow. The topographic variables of local slope and upslope contributing area, calculated by the model, can also be used to determine the discharge per unit width, q, or the average flow velocity, V , at any point in the catchment, thus defining the three-dimensional runoff processes affecting erosion and sediment transport. The simplest case that can be considered is that of steady-state runoff with a uniform rainfall excess throughout the catchment. These assumptions are adopted here for simplicity to illustrate the method. For this assumption q = rA/b, where r is the rainfall excess rate, A is the upslope contributingarea and b is the width of the contour element, asdefined in Figure 1 . The flow velocity can then be calculated from the estimates of q, local slope, and surface roughness using Manning's equation or a similar equation (Henderson, 1966) , as deemed appropriate. A full dynamic analysis is much more complex. An example of the potential relative erosion and deposition rates for the Geebung Creek catchment predicted by combining these concepts is presented in Figure 11 . In making these predictions it was assumed that sediment transport was transport capacity limited and not detachment limited. The effects of vegetation on erosion and deposition were also neglected. Because of these simplifications the results are presented in terms of potential relative erosion and deposition rates, rather than absolute values, and so only reflect the influence of topography on erosion and deposition. In Figure 11 erosion is negative and deposition is positive. Yang's (1973) unit stream power equation was chosen as the sediment transport capacity equation in this analysis because it is computationally simple, although any one of several equations that have been proposed in the literature (e.g. Alonso et al., 1981) could be used. This equation can be written as where C is the sediment concentration, P is the unit stream power, P,, is the critical unit stream power at incipient sediment motion, and y and p are constants that are functions of the median sediment size, the kinematic viscosity of water, and the terminal fall velocity of sediment particles in water. Unit stream power is defined as the time rate of potential energy dissipation per unit weight of water (= VS, where V is the flow velocity and S is the slope). Moore and Burch (1986a, 1986b) developed equations for estimating P for sheet and rill flow and demonstrated that equation (8) reliably predicts the sediment transport capacity of shallow sheet and rill flow for medium to coarse textured non-cohesive soils as well as finely aggregated clay soils (Moore and Burch, 1986a) . As an example, the equation for shallow sheet flow is where p is the density of water. Combining equations 8,9 and 10 allows the sediment flux per unit width to be calculated at any point in a catchment in terms of the soil properties P,,, n, y and fl (constants) and the hydrologic and topographic variables q and dylds, respectively, by making the appropriate assumptions (i.e. steady-state runoff). More detail is given by Moore and Burch (1986b) . As stated previously, most physically based erosion models developed to date represent erosion as a twodimensional process. However, Moore and Burch (1986b) showed that convergence and divergence of overland flow can significantly affect erosion and deposition. Therefore, the ability to realistically represent the effects of three-dimensional terrain on surface runoff, and hence erosion and deposition, using the outputs from the topographic model provides a significant improvement in our predictive capabilities. The inclusion of three-dimensionality into erosion models also allows such features as the transition of sheet flow to rill flow and the resulting geometry of rill networks to be better represented and predicted. The increase in sediment transport capacity as overland flow is concentrated in rill and channel networks and the influence of surface cover on soil properties and soil detachments are currently recognized as important processes to be represented in predictive models. Erosion and deposition models that realistically represent the three-dimensionality of natural terrain can be used to identify areas of high potential erosion that require special conservation treatment without having to resort to the treatment of entire catchments. This would allow for more efficient and effective use of the limited financial resources available for reducing erosion and sediment discharges from catchments. A model for determining the topographic attributes of slope, aspect, and upslope contributing area distributed across three-dimensional terrain is described. The discretization of the terrain produces irregularly shaped polygons based on contour lines and their orthogonals, streamlines, or flow lines. This form of discretization gives the model numerous advantages over existing models, particularly when the model is applied to problems involving water movement over and through natural three-dimensional terrain. The model was originally developed for hydrological applications and as a result its structure is well suited to representing natural earth surface processes, especially those relating to the hydrology, geomorphology, and sedimentology of landscapes and the physiology of native vegetation. The algorithms used in the model are described according to the two core computer programs that progressively analyse input data that are usually in the form of a digitized topographic map. Topographic, geomorphic, soil, vegetative, and hydrologic attributes of a small catchment in southeastern Australia are used to demonstrate three applications of the model. In the first application the topographic attributes calculated by the model are used to predict the natural occurrence of seepage or saturation zones within the catchment as an aid in the interpretation of its hydrological behaviour. Secondly, the computed topographic attributes are used to predict the distribution of potential daily solar radiation across the catchment to help interpret the distribution and species composition of tree associations or groupings previously surveyed within the catchment. Finally, the potential relative erosion and deposition rates within the catchment are estimated by combining a sediment transport function with procedures for simulating the hydrology of three-dimensional land surfaces using the calculated topographic attributes. These applications of the model are introduced to demonstrate its versatility and by no means exhaust the possibilities for other applications. Brief description of a comprehensive three-dimensional process-response model of landform development Estimatingsediment transport capacity in watershed modelling A three-dimensional simulation of slope forms Altitudinal distribution of several eucalypt species in relation to other Band ANSWERS: amodel for watershed planning A physically based variable contributing area model of basin hydrology', Hydrological Sciences Binns, D. 1984. 'Vegetation studies', in Eden Catchment Project Comparative hydrologic behaviour of forested and cleared catchments in southeastern Australia PTLOC-A FORTRAN subroutine for determining the position of a point relative to a closed boundary Quantifyingsource areas through land surface curvature and shape The varying source area of streamflowfrom upland basins Mathematical model and concept of equilibrium in connection with slope shear ratio 12 I I -1 220, 1226. environmental factors in southern New South Wales Research Paper No A summary of some surface fitting and contouring programs for noisy data Applied Hydrology Effect of wildfire and logging on the hydrology of small catchments near Eden Modelling erosion and deposition: topographic effects'. Transactions of Arneriian Society qf Hydrologic characteristics and modellingofa Distribution of slope steepness in the Palouse region of Washington The extraction of drainage networks from digital elevation data', Computer Vision Graphics O'Loughlin, E. M. 198 I. 'Saturation regions in catchments and their relations to soil and topographic properties' A mathematical model of soil erosion and deposition processes: I An algorithm for computing slope and aspect from elevations Bioclimatic factors in forest site potential'. Forest Ecology and Management VIEWIT Computation of seen areas, slope and aspect for land use planning'. Pacifrc Southwest Forest and Range Experiment Station Incipient motion and sediment transport Surface hydrology: 1. Explanation of phenomena, 11. Distribution ofraindrops, 111. Causes of lateral flow, IV. Flow in sloping layered soil, V. In-surface transient flow small forested catchment in southeastern New South Wales. Pre-logging condition The authors wish to thank Dr. Mike Hutchinson, CSIRO Division of Water and Land Resources for reviewing an early draft of this paper and for his many comments and suggestions during the course of the study. This work was partially supported by a grant (Project 85/131) from the Australian Water Research Advisory Council.