To be printed in Proceedings of the International Congress on Modelling and Simulation (MODSIM 2001). Evolution of a Hillslope Hydrologic Model F.G.R. Watson a, R.B. Grayson b, R.A. Vertessy c, M.C. Peel b and L.L. Piercea a Watershed Institute, Cal. State Univ. Monterey Bay, Seaside, CA, 93955, USA (fred_watson@csumb.edu) b CRC for Catchment Hydrology, Dept. of Civil and Environmental Engineering, Univ. of Melbourne, Parkville, VIC, 3052, Australia c CRC for Catchment Hydrology, CSIRO Land and Water, Canberra, ACT, 2601, Australia Abstract: The division of large landscapes into ‘hillslopes’ is an attractive concept in hydrologic modelling. It essentially reduces the dimensionality of the landscape from three dimensions to two: vertical and lateral. This results in economy with respect to both conceptual and computational complexity. The evolution of a hillslope hydrologic model is described, from beginnings in the traditional distribution function paradigm, to the presentation of a new ‘catena routing model’ (CRM). This development has been achieved adaptively. Each new application of the model to a different study area has revealed limitations in its predecessor, and has necessitated an ongoing search for a practical, general hillslope model. Keywords: Large scale; Hillslope; Catena; Hydrology; Distribution modelling; Macaque; Topmodel 1. INTRODUCTION Hydrologic models differ in the way lateral flow processes are represented, from not at all (i.e. lumped models), through semi-distributed models with lateral transfers at scales larger than hillslopes, to distribution-style conceptualisations (where the ‘elements’ are not spatially contiguous in the landscape, but rather are in some distribution based on process understanding), to fully explicit models. Representation of spatial heterogeneity within models is useful because it reduces the scale of representation of processes down to a level that is more commensurate with the spatial scale of our understanding of the processes themselves, and the dominant controls on these processes, such as land use change. In many cases, this scale is the ‘hillslope’ scale [Fan and Bras, 1995], that of loworder watersheds and the division of these into xeric through mesic areas. Lumped and semidistributed models are not suitable for investigations relating to processes operating at this scale. At the opposite end of the continuum, fully explicit models often represent more than is needed, and are data-demanding. Models that based their spatial structure around hillsopes, and the xeric-mesic catenae within them attempt to find the middle ground that is optimal for a certain type of landscape investigation. The Macaque model is one such model [Watson et al., 1998, 1999b; Peel et al., 2000a, 2001]. We describe the evolution of its hillslope sub-model. 2. HORIZONTAL PARTITIONING A central feature of Macaque is the parsimonious representation of spatial heterogeneity within large watersheds. For this, the model adopts the spatial structure set forth in the RHESSys model by Band et al. [1993]. Large watersheds are divided into many “hillslopes”, these being the segments of land between each node in a river network derived from analysis of a digital terrain model (DTM). This captures broader scales of heterogeneity, such as where different hillslopes are subject to differing solar radiation regimes. Within hillslopes, the major remaining heterogeneity is often aligned along a catena from ridge-top to valley-bottom. In reality, such catenae are often not ideal sequences, aligned along a single flow path down the hillslope. However, wetness indices can be calculated which tend to reflect the position of each part of the landscape along a hypothetical, continuous flow path. The TOPMODEL wetness index, ln(a/tanβ), [where a is upslope area per unit contour length and β is slope; Beven, 1997] is the most widely used for this purpose. Macaque divides each hillslope into about 10 different groups of cells, each associated with a different wetness index value. These groups are the elementary spatial units (ESUs) of the model. To be printed in Proceedings of the International Congress on Modelling and Simulation (MODSIM 2001). 3. LATERAL SUBSURFACE FLOW – THE TOPMODEL DFM The next step is to consider how to represent the lateral movement of water between ESUs. Three schemes have been used within Macaque to date. The model was originally designed for use in forest situations with soils of very high infiltration capacity. So overland flow down hillslopes was ignored. Most lateral redistribution of water in these conditions occurs as subsurface flow. The first scheme implemented was what we refer to as the distribution function modelling (DFM) approach, as described by Beven [1997] for single hillslopes and by Band et al. [1993] and others for large landscapes comprised of multiple hillslopes. This assumes a steady state flow through hillslopes such that the ‘shape’ of the water table when viewed as a section down the catena remains constant for a given rate of recharge, and varies as the entire hillslope becomes wetter through a higher rate of recharge, and consequently, a higher total saturated soil moisture storage. The approach assumes a spatially constant rate of recharge, r, and a spatially constant vertical saturated soil hydraulic conductivity (K) profile. If simple mathematical formulae are used to represent the K profile, then the depth to the water table, zx, of any point, x = [0…L], in the landscape can be predicted by an analytical mathematical solution to the governing equations. For example, the traditional assumption is for an exponential decay in Kx: K x (z x ) = e − f zx (1) where f is a shape parameter. Equation (1) can be expressed in terms of soil saturation deficit, Sx: K x (S x ) = e − S x m (2) where m is also a shape parameter that is related to f by way of the soil active porosity, ∆θ,x, which is assumed to be vertically constant: m= ∆θ , x , f S x = z x ∆θ , x (3) This leads to the following distribution function for the local water table depth Sdist,x [Beven and Kirkby, 1979]: ( S dist , x = S − m wx − w wx = ln (a x tan β x ) ) (4) (5) where wx is a wetness index at point x, ax is upslope area per unit contour length, and βx is terrain slope. Models using steady-state expressions of this form are termed distribution function models (DFMs). Their prediction of the lateral distribution of moisture deficits requires only knowledge of landscape attributes (ax and βx), and accounting of the mean hillslope moisture deficit, S . Models of this form instantly redistribute water each time step to maintain the shape of the deficit curve defined by the distribution function. This can be thought of as an “infinite lateral conductivity” assumption – i.e. there is no time delay in lateral redistribution. Further, Beven and Kirkby’s analytical framework conveniently yields a term for the subsurface (baseflow) discharge of the entire hillslope, qbase, based only on the hillslope moisture deficit and a scaling parameter q0: ( qbase = q0 exp − S m 4. ) (6) AN EXPLICIT DISTRIBUTION FUNCTION MODEL (EDFM) For practical applications, where we are interested in being able to test the location of runoff producing areas in the field, we developed a more explicit model of runoff generation [Watson et al., 1998]. Under this ‘explicit DFM’ (EDFM), baseflow, qbase,i, for each ESU calculated separately as follows. Firstly, the ESU is thought of as having two horizontal portions: the ‘saturated portion’, pi : 0 ≤ pi ≤ 1, that is saturated to the surface, and the ‘unsaturated portion’, (1 pi), whose water table lies somewhere below the surface. These portions are assumed to vary with the saturation deficit of the ESU as a whole (the saturation deficit represents the location of the water table, which divides the ESU vertically into saturated and unsaturated zones). By assuming a uniform surface gradient, ∆e,i/∆x, underlain by a uniform water table gradient, we can estimate the saturated portion as a linear function of saturation deficit with the help of a parameter, Sq,i, that specifies the saturation deficit at which the saturated portion is exactly zero: pi = S q ,i − Si (7) S q ,i Baseflow for each ESU is assumed to exfiltrate through the saturated portion of the ESU at a rate controlled by the surface K, Ksurf,i, and the hydraulic gradient, ∆h,i/∆x, of the ESU: qbase ,i = pi K surf ,i ∆ h ,i ∆x (8) To be printed in Proceedings of the International Congress on Modelling and Simulation (MODSIM 2001). where the hydraulic gradient is estimated as the surface gradient, ∆e,i/∆x, scaled by a parameter, ρ: 0 ≤ ρ ≤ 1: ∆ h ,i ∆ = ρ e ,i ∆x ∆x (9) For generalized DFMs, where the gradient of saturation deficit with respect to wetness index is specified as a key parameter, ∆S/∆w, (see Equation (18) below), Sq,i can be estimated from ∆S,i, the saturation deficit at the dry end of an ESU when the wet end is just becoming saturated, and ∆w,i, the range of wi values observed across all cells within an ESU [Watson, 1999]: S q ,i = ∆ S ,i 2 = − 1 ∆S ∆ w, i 2 ∆w (10) In catena routing models (see below), a more certain estimation is possible: S q ,i = ∆θ ,i (1 − ρ ) ∆ e ,i 2 Firstly, in analytical solutions, the wetness index is governed by the assumed shape of the Kx(zx) function. Ambroise et al. [1996] give the wetness indices and distribution functions that result from some alternate assumptions about the K profile. It is important to realise that in order for these analytical solutions to exist, all expressions of K versus depth must have an inversion that can be integrated. This is restrictive. Observed profiles of K [e.g. Davis et al., 1999] are often best represented by more complex equations representing the sum of two or more simple equations. For example, we have found the following general form to be able to characterize a variety of (in particular) deep soils: Ki (11) The saturated proportion is also used in the estimation of stormflow runoff generated by rain falling directly on saturated areas. Maps can thus be drawn of the total runoff (qΣ) producing areas for each time step of simulation (by colouring each ESU according to its runoff), and of the surfce saturation of the hillslope (by colouring each ESU according to its saturated portion). We prefer this scheme because it forces a reconciliation between the dynamics of predicted water table shapes, predicted spatial distribution of runoff producing areas, and predicted hydrographs. We deem the model valid only if all three are satisfactorily predicted. In the Ettercon3 catchment, this was determined using piezometer nests, field survey of source areas, and inference of the existence of variable source areas from the dependence of the storm flow / precipitation ratio on antecedent conditions [Campbell, 1998; Watson, 1999]. 5. which these can be manipulated. As explained below, the wetness index is restricted by the need for the inversion to be integrable, and the distribution function is restricted by need to maintain continuity. GENERALIZING THE DISTRIBUTION FUNCTION MODEL (GDFM) Early experiments with Macaque [Watson et al., 1996] showed that by varying the TOPMODEL ‘m’ parameter, the model was able to reproduce either the observed hydrograph or the observed catena profile of water table depth [Campbell, 1998], but not both for the same value of m. Initially, this caused us to question the form of the wetness index and the distribution function. However, there are limitations to the degree to  K min,i + (K surf ,i − K min,i )e − fz = 0  zi < zmax,i otherwise (12) This satisfies the observed tendency for conductivity to dramatically decline in soil’s Chorizon, as well as the tendency for very high K values in the O and A horizons of many forest soils [Davis et al., 1999]. This equation cannot be inverted to produce an analytically integrable equation for z in terms of K (or, analogously, S in terms of K). Thus, a wetness index and distribution function cannot be calculated from this equation in the same way as it could for the K equations presented by Ambroise et al. [1996]. This represents a limitation in the range of wetness indices available to us that have matching K profiles. Secondly, given a particular wetness index, the way in which we can use this index in a distribution function is also limited. It must provide continuity of water balance. The essential requirements of a distribution function are that it predicts local moisture conditions based on mean moisture conditions and a measure of location (i.e. a wetness index): ( Si = f S , wi ) (13) Such an equation allows us to specify local moisture conditions based on global moisture accounting. However, if we are to do this, we must ensure that the sum of local moisture values is the global value: To be printed in Proceedings of the International Congress on Modelling and Simulation (MODSIM 2001). [Campbell, 1998]. This cannot be represented by a GDFM. 1 n ∑ [Si Ai ] = S ΣA i =0 (14) where Ai is the area of ESU i, and ΣA is the total area of the hillslope. The only distribution function that satisfies this is one of the form: ( ) f S , wi = S + g (wi ) (15) where g() is some function of wetness index alone [Watson, 1999]. The required equality is then: 1 n ∑ [ g (wi )Ai ΣA i = 0 ] = 0 (16) The most general form of g that satisfies Equation (16) is: g (wi ) = ( ∆S wi − w ∆w ) (17) where ∆S /∆w is the gradient of saturation deficit with respect to wetness index (negative). Thus, the generalised distribution function (GDF) is thus: S dist ,i = S + ( ∆S wi − w ∆w ) (18) This analysis thus shows that the original Topmodel distribution function (Equation (4)) is a member of a rather exclusive set of functions that can be used to achieve steady state distribution function modelling. Given the requirements of matching K functions, and continuity of water balance, we are only free to choose from wetness indices derived from invertible K functions, and different values of a scaling parameter ∆S/∆w. Interpreting this in physical terms, a GDFM (e.g. Topmodel) can only simulate profiles of moisture deficit down a hillslope catena that are of fixed shape, moving vertically only in unison. Any increase or decrease in total hillslope moisture must involve exactly the same increase or decrease in moisture at all parts within the hillslope. This was expressed by Beven and Kirkby as the ‘uniform recharge assumption’. In shallow soils with uniform vegetation cover, such an assumption may hold. But in deep forest soils with welldefined riparian zones, non-uniform recharge is ubiquitous. Data from the Maroondah catchments show that deep (>20 m) water tables away from the stream move very little over time, while shallow water tables (<2 m) adjacent to saturation areas move both seasonally and with every storm 6. LIMITING THE DISTRIBUTION FUNCTION MODEL (LDFM) The above discussion is underlain by the paradigm that the only means by which we allow local saturated moisture deficits to vary is by way of the distribution function. In order to escape this paradigm, a means is required by which local moisture accounting can be combined with distribution function accounting in such a way that continuity of water balance is not violated. Watson et al. [1998] implemented a solution termed ‘limited distribution function modelling’ (LDFM). Under this method each ESU accepts a lateral flux to or from the remainder of the hillslope by way of a distribution function that is limited by an additional governing parameter. Under the conventional GDF paradigm, the net lateral flux from ESU i at time t, qlaterai,i,t, would be: qlateral ,i ,t = S dist ,i ,t − Si ,t −1 (19) The LDFM approach adds an additional factor, δ: 0 ≤ δ ≤ 1: qlateral , i , t = (S dist , i , t − Si , t −1 ) δ (20) This effectively removes the steady state assumption and slows down the imposition of temporally averaged moisture patterns imposed by the distribution function, allowing local water tables to fluctuate about this average with each storm. Because under the EDFM approach, streamflow is directly tied to surface saturation, we have found that by calibrating this factor, we can simulate realistic water table dynamics for both shallow and deep water tables, concurrent with accurate streamflow dynamics. Both of these are crucial, as we are interested in the control of tree water consumption by soil moisture dynamics, and in turn, the effect of tree water consumption on the flow duration curve and risks to domestic water supply [Watson et al., 1999b]. Because the limiting factor, δ, is applied uniformly across the hillslope to all lateral fluxes imposed by the distribution function, continuity of water balance is assured. We have found that the optimal values of δ calibrated against both flow and piezometer data can be very low, in the order of 0.001 [Watson, 1999]. This implies that the “infinite lateral To be printed in Proceedings of the International Congress on Modelling and Simulation (MODSIM 2001). conductivity assumption” does not hold and that in fact lateral redistribution is very slow. It also implies that the exact choice of wetness index is not particularly critical, so long as it predicts longterm spatial moisture patterns. Further, since local moisture, and hence streamflow is now controlled more by local vertical fluxes, there is a diminished need to require that the K profile equation and the distribution function can be analytically linked (as in Beven and Kirkby, 1979). Thus the two main limitations of the purely GDFM approach are somewhat relaxed by the adoption of the LDFM approach. Yet the parsimony and computational efficiency of distribution function modelling in general is retained, with the addition of a single parameter, δ. 7. SPATIAL PARAMETERISATION AND THE REALITY OF APPLIED MODELLING Landscapes are not managed at the hillslope scale. Rather, they are managed at a regional level, as collections of a large number of hillslopes. Thus while the hillslope is a convenient unit for understanding processes such as runoff generation, we must simulate regions in order to understand the interplay between hillslope processes and spatially varying climate, vegetation, soils, and management policy. The parameters of a hillslope model are dependent upon such factors, and so a ‘regionalisation’ scheme must be developed to estimate values of these parameters across a large spatial domain. Topmodel and other GDFMs have two degrees of freedom, an equation for wetness index, w, and a scaling parameter (either m or f). EDFMs have an additional parameter ρ, and LDFMs have a further parameter, δ. None of these parameters can be measured effectively in the field. They must be calibrated against hydrographic and water table data, or estimated through regionalisation of previously calibrated values. Watson et al. [1998, 1999b] calibrated ρ and δ against both flow and pizeometer data in the Maroondah catchments. Peel et al. [2001] and [Peel, 1999] calibrated ρ against flow data in a number of diverse catchments throughout Australia. Returning to both Maroondah and the Thomson catchments, Peel et al. [2000a,b] again calibrated ρ against gauging data, and also developed a simple scaling relationship between ρ and catchment size. The Maroondah soils are relatively homogeneous, and many field data on runoff processes are available, so we are comfortable that calibrated values of ρ and δ are robust when applied to ungauged subcatchments. However, when soil-type varies, as is often the case within large landscapes, values of ρ and δ become uncertain and the approach falters. Further, when variation in soil parameters occurs within a hillslope, distribution function models in general are of limited use. 8. ABSTRACTION OF A GEOMETRIC CATENA FROM A HILLSLOPE In searching for the optimal, practical hillslope model, we revist their initial utility: the parsimony achieved through representation of landscape spatial structure as a collection of catenae through which flow moves in two dimensions (vertical and downslope), rather than as a three dimensional network. The catena is a useful abstraction of a hillslope, encapsulating both the boundary conditions and dominant flow directions of a hillslope. This abstraction is the primary benefit of hillslope/catena modelling. Simulation of lateral redistribution through the use of distribution functions is secondary, and not essential. An alternative is to simulate lateral flow explicitly. This is the approach currently adopted within Macaque. Every time step, downslope flux is calculated for each point along the catena using Darcy’s Law. Water table levels respond to these flows, and surface flow is produced when the water table rises to the surface. The required parameters are all obtainable from terrain and soils data. The critical requirement is a representation of catena geometry suitable for explicit subsurface flow calculations. A method is proposed here for constructing a catena geometry that is representative of the important geometry of hillslopes characterised by gridded terrain data. Firstly, we use a wetness index to locate all hillslope cells along a moisture gradient that will be used as a proxy for the soil catena. The gradient is discretised to isolate groups of cells as distinct ESUs. The catena geometry is constructed as a series of trapezoidal prisms, each a direct analogue of a corresponding ESU along the wetness gradient. Each trapezoid is parameterised by its area, Ai, its width at the lower boundary, Wi, and its convergence, ci. Convergence is here defined as the ratio of upper to lower boundary widths of the trapezoidal prisms: ci = Wi−1 Wi (21) To be printed in Proceedings of the International Congress on Modelling and Simulation (MODSIM 2001). The trapezoid length, Li, can be derived from these parameters by a re-arrangement of the formula for the trapezoid area: A 2 Li = i Wi 1 + ci (22) The required parameters, Ai, Wi, and ci are calculated through terrain analysis of a gridded digital elevation model (DEM). Ai is simply the sum of the area of all cells in each ESU. Width calculations are propagated up from the ESU at the bottom of the hillslope using Equation (21). The lower boundary width, Wi, of the lowest ESU can be estimated from the widths and number of stream cells below the hillslope, which are in turn defined as cells with an upslope area greater than some threshold determined from topographic maps or aerial photos. The convergence, cj, of a single cell, j, is estimated from its valency, vj, which is the number of upslope neighbour cells for which cell j is their steepest descent neighbour, with the exception that cells with no upslope neighbour have a convergence of one: c j = max (1, v j ) (23) Mean ESU convergence, ci, is calculated as the mean convergence of all cells in an ESU. 9. A CATENA ROUTING MODEL (CRM) Armed with the above parameters, and techniques for estimating them from terrain data, we are able to devise a simple catena routing model (CRM). Each time step, the lateral subsurface flow leaving each ESU (per unit area of the ESU) is calculated using Darcy’s Law: qlateral ,i = Ti ( zi ) ∆ h ,i Li Li (24) where Ti() is a transmissivity function specific to each ESU, evaluated at water table depth, zi, calculated by integrating Equation (12), T (z ) = z max ∫ K ( z ) dz z =0 (25) and ∆h,i is the hydraulic head difference between the ESU and its downslope neighbour: ∆ h,i = ∆ e,i + ∆ z ,i (26) ∆e,i is an estimate of elevation difference between the two ESUs, and is calculated from the surface gradient and length of the ESU: ∆ e ,i = Li tan (β i ) (27) where βi is the mean surface slope of the ESU, and the ESU length, Li, is derived from Equation (22) The actual difference in mean elevations for the two ESUs cannot be used, because the ESU that is the drier of the two with respect to the wetness index gradient is not guaranteed to be higher in elevation. Similarly, the elevation range amongst all cells within an ESU cannot be used, as this can be zero. Therefore we use the effective hydraulic gradients across an ESU to represent that between the ESU and its neighbour. ∆z,i is the difference in water table depths, zi, between the two ESUs, with zi calculated from Equation (3). During development of the CRM algorithm, a simplified version was used in Macaque applications to the Thomson catchment in Australia [Peel et al., 2000a,b], to eight diverse catchments from around Australia [Peel 1999; Peel et al., 2001], and to the Salinas Valley in California [Watson et al., 1999a]. The simplification assumed a square hillslope with uniform width, W, estimated as: Wi = W = ΣA (28) and the trapezoids were in fact rectangles with length: Li = Ai W (29) While these applications were not able to represent hillslope convergence, they are appealing because they test the approach’s ability to respond directly to measurable soil parameters. For example, in the Salinas application, regional soil maps were used to represent highly variable soil depths within the landscape, and within hillslopes themselves. The model was thus able to simulate observed phenomena such as depressed water tables beneath sandy alluvium at the base of hillslopes with otherwise skeletal soils. The full convergence algorithm, using Equation (22) for trapezoid length, is currently being tested. 10. CONCLUSIONS Pure distribution function models (DFMs) such as Topmodel are theoretically attractive, but in many situations do not fully reproduce observed spatial runoff dynamics due to restrictions imposed by their basic structure. Explicit and limited DFMs can represent a wider range of dynamics, but To be printed in Proceedings of the International Congress on Modelling and Simulation (MODSIM 2001). require additional, sometimes uncertain parameters to do so. A catena routing model (CRM) is outlined that attempts to solve these problems by introducing a novel method of abstracting catena geometry from cell-based hillslope data, and using more easily defined parameters to support a Darcian flow algorithm within the abstract catena. 11. ACKNOWLEDGEMENTS For fertile discussion over the years, we thank Rob Lamb, Larry Band, Scott Mackay, Richard Lammers, Christina Tague, Keith Beven, and Andrew Western. 12. REFERENCES Ambroise, B., K. Beven, and J. Freer, Toward a generalisation of the TOPMODEL concepts: Topographic indices of hydrological similarity. Water Resour. Res., 32:7:2135-2145, 1996. Band, L.E., P. Patterson, R. Nemani, and S. Running, Forest ecosystem processes at the watershed scale: incorporating hillslope hydrology. Agric. For. Meteorol., 63:93-126, 1993. Beven, K.J. TOPMODEL: A critique. Hydrol. Proc., 11:1069-1085, 1997. Beven, K.J. and M.J. Kirkby, A physically based, variable contributing model of basin hydrology. Hydrol. Sci. Bull., 24:43-69, 1979. Campbell, R., Groundwater Behaviour in a Small Forested Catchment. MSc, Hydro. and Env. Sci., U. Melbourne, 1998. Davis, S.H., R.A. Vertessy, and R. Silberstein, The sensitivity of a catchment model to soil hydraulic properties obtained by using different measurement techniques. Hydrol. Proc., 13:677688, 1999. Fan, Y. and R.L. Bras, On the concept of a representative elementary area in catchment runoff. Hydrol. Proc., 9, 821-832, 1995. Peel, M.C., Annual runoff variability in a global context. PhD, Geog. and Env. Stud., U. Melbourne, 1999. Peel, M.C., T.A. McMahon, B.L. Finlayson, and F.G.R. Watson, Identification and explanation of continental differences in the variability of annual runoff. J. Hydrol., 250:224-240, 2001. Peel, M.C., F.G.R. Watson, and R.A. Vertessy, The impact of logging and fire on water yield: Predicting the Thomson Catchment by Macaque. Water, 27:20-24, 2000a. Peel, M.C., F.G.R. Watson, R.A. Vertessy, J.A. Lau, I.S. Watson, M.W. Sutton, and B.G. Rhodes, Predicting the water yield impacts of forest disturbance in the Maroondah and Thomson catchments using the Macaque model. View publication stats Rep. 00/14, CRC for Catch. Hydrol., Vic., Australia, 2000b. Watson, F.G.R., Large scale, long term, physically based modelling of the effects of land cover change on forest water yield. PhD, Civil and Env. Eng., U. Melbourne, 443 pp. earthsystems.csumb.edu/~fwatson, 1999. Watson, F.G.R., R.B. Grayson, R.A. Vertessy, and T.A. McMahon, Large scale distribution modelling and the utility of detailed ground data. Hydrol. Proc., 12:873-888, 1998. Watson, F.G.R., L.L. Pierce, M. Mulitsch, W. Newman, J. Nelson, and A. Rocha, Spatial modelling of the impacts of 150 years of land use change on the carbon, nitrogen, and water budgets of a large watershed. ESA Annual Meeting, 1999, Spokane, USA, 1999a. Watson, F.G.R., R.A. Vertessy, and L.E. Band, Distributed parameterization of a large scale water balance model for an Australian forested region. IAHS Pub. 235:157-166, 1996. Watson, F.G.R., R.A. Vertessy, and R.B. Grayson, Large scale modelling of forest hydrological processes and their long term effect on water yield. Hydrol. Proc., 13:689-700, 1999b.