r JOURNAL VOL. 34, NO. 1 OF THE AMERICAN WATER RESOURCES ASSOCIATION AMERICANWATERRESOURCES ASSOCIATION FEBRUARY 1998 A LARGE AREA HYDROLOGIC PART I: MODEL MODELING AND ASSESSMENT DEVELOPMENTl J. G. Arnold, R. Srinivasan, R. S. Muttiah, and J. R. Williams2 ABSTRACT: A conceptual, continuous time model called SWAT (Soil and Water Assessment Tool) was developed to assist water resource managers in assessing the impact of management on water supplies and nonpoint source pollution in watersheds and large river basins. The model is currently being utilized in several large area projects by EPA, NOM, NRCS and others to estimate the off-site impacts of climate and management on water use, nonpoint source loadings, and pesticide contamination. Model development, operation, limitations, and assumptions are discussed and components of the model are described. In Part II, a GIS input/output interface is presented along with model validation on three basins within the Upper Trinity basin in Texas. (KEY TERMS: simulation; surface water hydrology; erosion; sedimentation; nonpoint source pollution; large area modeling; plant growth; agricultural land management.) reasonable results. The model must correctly reflect changes in land use and agricultural management on stream flow and sediment yield. Available models with these capabilities are generally limited by spatial scale: Available river-basin models generally do not link outputs to land use and management adequately to evaluate management strategies. Also, most are single-event models. We chose good agricultural management models to link with simple, efficient, yet realistic routing components for the purpose of capturing management effects on large river basins through long-term simulations. The objective of this overview is to briefly describe an overview of model operation, model applications, and a description of model components of a river basin scale model called SWAT (Soil and Water Assessment Tool). INTRODUCTION Large area water resources development and management require an understanding of basic hydrologic processes and simulation capabilities at the river basin scale. Current concerns that are motivating the development of large area hydrologic modeling include climate change, management of water supplies in arid regions, large scale flooding, and off site impacts of land management. Recent advances in computer hardware and software including increased speed and storage, advanced software debugging tools, and GIS/spatial analysis software have allowed large area simulation to become feasible. The challenge then is to develop a basin-scale model that: (1) is computation ally efficient; (2) allows considerable spatial detail; (3) requires readily available inputs; (4) is continuous-time; (5) is capable of simulating land-management scenarios; and (6) gives LITERATURE REVIEW Integrated water management of large areas should be accomplished within a spatial unit (the watershed) through modeling. Integrated water management can be viewed as a three or more dime~sional process centered around the need for water, the policy to meet the needs, and the management to implement the policy. Watershed modeling is fundamental to integrated management. Watershed models abound in the hydrological literature (Singh, 1989) and state-of-the-art of watershed modeling is reasonably advanced. However, these models have yet to become common planning or decision making tools. A lPaper No.96089 of the Journal of the American Water ResourcesAssociation. Discussions are open until October I, 1998. 2Arnold and Williams, Agricultural Engineer and Hydraulic Engineer, respectively, USDA.Agricultural Research Service, 808 East Blackland Road, Temple, Texas 76502; Srinivasan and Muttiah, Associate Research Scientists, Texas Agricultural Experiment Station, 808 East Blackland Road, Temple, Texas 76502 (e.m (Arnold): arnold@brcsuncp.tamu.edu). JOURNAl OF THE AMERICAN WATER RESOURCES ASSOCIATION 73 JAWRA 1t Arnold, Srinivasan, Muttiah, and Williams majority of watershed models simulate watershed response without or with inadequate consideration of water quality. If these models are to be used for environmental or ecological modeling, they must consider water quality (Singh, 1995). After the development of the Stanford Watershed Model (Crawford and Linsley, 1966) numerous operational, lumped or "conceptual" models have been developed. These include: SSARR (Rockwood et al., 1972), the Sacramento model (Bumash et al., 1973), the tank model (Sugawara et al., 1976); HEC-1 (Hydrologic Engineering Center, 1981), HYMO (Williams and Hann, 1973), and RORB (Laurenson and Mein, 1983). In these models, some processes are described by differential equations based on simplitied hydraulic laws, and other processes are expressed by empirical algebraic equations. More recent conceptual models have incorporated soil moisture replenishment, depletion and redistribution for the dynamic variation in areas contributing to direct runoff. S~veral models have been developed fro II! this concept which use a probability distribution of soil moisture including the ARNO model (Todini, 1996; Zhao, 1984; Moore and Clarke, 1981) or the use of a topographic index, as in TOPMODEL (Beven and Kirk by, 1979; Beven et al., 1984). Jayatilaka et. al. (1996) recently developed a variable source conceptual model that shows promise for incorporation into comprehensive models. Another class of hydrological models is a differential model based on conservation of mass, energy, and momentum. Examples of differential models include SHE (Abbott et al., 1986a, 1986b), IDHM (Beven et al., 1987), and Binley et al. (1989). The SHE model simulates water movement in a basin with the finite difference solution of the partial differential equations describing the processes of overload and channel flow, unsaturated and saturated subsurface flow, interception, ET, and snowmelt. The spatial distribution of catchment parameters is achieved by representing the basin on an orthogonal grid network. Jain et al. (1992) successfully applied the SHE model to an 820 km2 catchment in central India. However, they note that the data requirements are substantial. Jain et al. (1992) also concluded that the strength of differential models like SHE "lies beyond the field of pure rainfall-runoff modeling, for which purpose traditional and simpler hydrologic models often perform equally well." In the early 1970s work also began on non-point source modeling in response to the Clean Water Act. The CREAMS model (Knisel, 1980) was developed to simulate the impact of land management on water, sediment, nutrients, and pesticides leaving the edge of a field. Several field scale models evolved from the JAWRA original CREAMS to simulate pesticide ground water loadings (GLEAMS; Leonard et at., 1987) and to simulate the impact of erosion on crop production (EPIC; Williams et at., 1984). Other efforts evolved to simulate hydrology and water quality of complex watersheds with varying soils, land use, and management. Several models were developed to simulate single storm events using a square grid representation of spatial variability (Young et at., 1987; Beasley et at., 1980). These models did not consider subsurface flow, ET or plant growth. Continuous models were also developed (Johansen et at., 1984; Arnold et at., 1990) but generally lacked sufficient spatial detail. SCALING ISSUES MODEL OPERATION SWAT is an operational or conceptual model that operates on a daily time step. The objective in model development was to predict the impact of management on water, sediment and agricultural chemical yields in large ungaged basins. To satisfy the objective, the model (a) does not require calibration {calibration is not possible on ungaged basins); {b) uses readily available inputs for large areas; {c) is computation ally efficient to operate on large basins in 74 JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATK)N Large Area Hydrologic Modeling and Assessment -Part a reasonable time, and (d) is continuous time and capable of simulating long periods for computing the effects of management changes. A command structure is used for routing runoff and chemicals through a watershed similar to the structure ofHYMO (Williams and Hann, 1973). Commands are included for routing flows through streams and reservoirs, adding flows, and inputting measured data on point sources (Figure 1). Using the routing command language, the model can simulate a basin subdivided into grid cells or subwatersheds. Additional commands have been developed to allow measured and point source data to be input to the model and routed with simulated flows. Although the model operates on a daily time step and is efficient enough to run for many years, it is intended as a long term yield model and is not capable of detailed, single-event flood routing. I: Model Development Q = (R- 0.28)2 R + 0.8 8 : R > 0.28 (2) R ~ 0.2 s Q=O.O, where Q is the daily surfa~e runoff (mm), R is the daily rainfall (mm), and s is a retention parameter. The retention parameter, s, varies (a) among watersheds because soils, land use, management, and slope all vary and (b) with time because of changes in soil water content. The parameter s is related to curve number (CN) by the SCS equation (USDA-SCS, 1972). 100 B = 254 (3) -1 CN The constant, 254, in Equation (3) gives s in mm. Fluctuations in soil water content cause the retention parameter to change according to the equation MODEL COMPONENTS 8= 81 The subbasin components can be placed into eight major divisions -hydrology, weather, sedimentation, soil temperature, crop growth, nutrients, pesticides, and agricultural management. 1- FFC FFC + exp[ wl -W2 ( FFC) ] (4) where SI is the value of s associated with CN 1, FFC is the fraction of field capacity, and wl and w2 are shape parameters. FCC is computed with the equation Hydrology FFC= The hydrology model is based on the water balance equation (Figure 2) swt = sw+ t L (Rj -Qi -ETi -Pi -QR;) (1) where SW is the soil water content minus the 15-bar water content, t is time in days, and R, Q, ET, P, and QR are the daily amounts of precipitation, runoff, evapotranspiration, percolation, and return flow; all units are in mm. Since the model maintains a continuous water balance, complex basins are subdivided to reflect differences in ET for various crops, soils, etc. Thus, runoff, is predicted separately for each subarea (Figure 3) and routed to obtain the total runoff for the basins. This increases accuracy and gives a better physical description of the water balance. Percolation. The percolation component uses a storage routing technique combined with a crack-flow model to predict flow through each soil layer. Once water percolates below the rootzone, it is lost from the watershed (becomes ground water or appears as return flow in downstream basins). The storage routing technique is based on the equation Surface Runoff Volume. Surface runoff is predicted for daily rainfall by using the SCS curve number ~quation (USDA-SCS, 1972) OF THE AMERICAN WATER RESOURCES ASSOCIATION (5) where SW is the soil water content in the root zone (mm) , WP is the wilting point water content(mm), (1,500 kPa for many soils), and FC is the field capacity water content (mm) (33 kPa for many soils). Values for wl and W2 are obtained from a simultaneous solution of Equation (4) according to the assumptions that s = S2when FCC = 0.6 and s = S3' when (SW-FC)/(POFC) = 0.5 There are two options for estimating the peak runoff rate -the modified Rational formula and the SCS TR-55 method (USDA-SCS, 1986). A stochastic element is included in the Rational equation to allow realistic simulation of peak runoff rates, given only daily rainfall and monthly rainfall intensity. t=l JOURNAL SW-WP FC-WP 75 JAWRA Arnold, Srinivasan, Muttiah, and Williams Figure 1. SWAT Model Operation Flow Chart. JAWRA 76 JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION Large Area Hydrologic Modeling and Assessment -Part I: Model Development Irrigation Typical Depths t.. Evaporation ~~ Layers OSOil ~/ Soil HOOt Soil Moisture &,,~ ~C'& ~ ~ ,".::==::::- -" Redistribution -:===:: .'$'v~ 0,," Percolate Profile ~o~ / Recharge ~ Revap 25m Return Flow ~ Shallow Transmission Losses Percolation from Shallow / Recharge to Deep Aquifer Aquifer Deep Aquifer Figure 2. Components of the Hydrologic Balance Simulated Within a SWAT Subbasin. where Hi is the hydraulic conductivity in mmh-l and FC is the field capacity minus wilting point water content for layer i in mm. The hydraulic conductivity is varied from the saturated conductivity value at saturation to near zero at field capacity. (6) where SWo and SW and the soil water contents (mm) at the beginning and end of the day, respectively; Llt is the time interval (24 h); and Tr is the travel time (h) through layer i. Thus, the percolation can be computed by subtracting SW from SW0. (9) where SCi is the saturated conductivity for layer i (mmh-l) and /3 is a parameter that causes Hi to approach zero as SW i approaches FCiThe equation for estimating /3 is (7) where O is the percolation rate (mm.d-l). The travel time, Tri, is computed for each soil layer with the linear storage equation -2.655 ~i = loglo TTi = (SWi -FCi) H, JOURNAL ("UL;" FCi ) (10) (8) I OF THE AMERICAN . The constant (-2.655) in Equation (10) was set to assure Hi = 0.002SCi at field capacity. Upward flow WATER RESOURCES ASSOCIATION 77 JAWRA Arnold, Srinivasan, Muttiah, and Williams Route Through Next Reach or R-rvoir Figure 3. Hydrologic Flow Chart of SWAT Subbasin Model. JAWRA 78 JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION Large Area Hydrologic Modeling and Assessment -Part may occur when a lower layer exceeds field capacity. Movement from a lower layer to an adjoining upper layer is regulated by the soil water to field capacity ratios of the two layers. Percolation is also affected by soil temperature. If the temperature in a particular layer is 0°C or below, no percolation is allowed from that layer. Lateral Subsurface Flow. Lateral subsurface flow in the soil profile (0-2 m) is calculated simultaneously with percolation. A kinematic storage model (Sloan et aZ., 1983) is used to predict lateral flow in each soil layer. q,at = 0.024 (2 S SC sin(a») ed L -rf -percgw -WUSA (13) Re O.8~a (14) (lO-e-Mt) -ea)1 )1 ra]] ra (15) The latent heat of vaporization, saturation vapor pressure, and slope of the saturation vapor pressure curve are all estimated with the temperature function (Arnold et al., 1993) The Hargreaves and Samani (1985) method estimates potential evapotranspiration as a function of extraterrestrial radiation and air temperature. Hargreaves' method was modified for use in SWAT by increasing the temperature difference exponent from 0.5 to 0.6. Also, extraterrestrial radiation is replaced by RAM X (maximum possible solar radiation at the earth's surface) and the coefficient is adjusted from 0.0023 to 0.0032 for proper conversion. The modified equation is where a is the constant of proportionality or the reaction factor. The relationship for water table height is (Arnold et aZ., 1993) e-aM + Evapotranspiration. The model offers three options for estimating potential ET -Hargreaves (Hargreaves and Samani, 1985), Priestley-Taylor (Priestley and Taylor, 1972~ and Penman-Monteith (Monteith, 1965). The Penman-Monteith method requires solar radiation, air temperature, wind speed, and relative humidity as input. If wind speed, relative humidity, and solar radiation data are not available (daily values can be generated from average monthly values), the Hargreaves or Priestley-Taylor methods provide options that give realistic results in most cases. (16) where V sa is the shallow aquifer storage (mm), Rc is recharge (percolate from the bottom of the soil profile) (mm), revap is root uptake from the shallow aquifer (mm), rf is the return flow (mm), percgw is the percolate to the deep aquifer (mm), WUSA is the water use (withdrawal) from the shallow aquifer (mm), and i is the day. Return flow from the shallow aquifer to the stream is estimated with the equation (Arnold et at., 1993): hi = hi-l (m above stream where Eo is evaporation (g m-2s-l), HV is latent heat of vaporization (J g-l ), ho is net radiation (J m-2s-l), O is slope of the saturation vapor density function (g m-3C-l ), S is soil heat flux (J m-2s-l), 'Yis psychometric constant (g m-3C-l), Pa is air density (g m-3), Cp is specific heat of air (J g-lC-l ), es is saturation vapor density (g m-3), ea is air vapor density (g m-3), r a is aerodynamic resistance for heat and vapor transfer (s m-l), and rc is canopy resistance for vapor transfer (s m-l). The Priestley-Taylor (1972) method provides estimated of potential evaporation based only on temperature and radiation is Ground Water Flow. Ground water flow contribution to total streamflow is simulated by creating a shallow aquifer storage. The water balance for the shallow aquifer is rfi = rfi e-at1t + Rc(1.0 -e-a,:\t) where h is the water table height, bottom), and J.Lis the specific yield. S)+PaCp(e. Eo = o(ho+ HV[o+y[(ra-rc where qlat is lateral flow (mm d-l), S is drainable volume of soil water (mh-l), <Xis slope (mm-l), ed is drainable porosity (mm-l ), and L is flow length (m). If the saturated zone rises above the soil layer, water is allowed to flow to the layer above (back to the surface for the upper soil layer). To account for multiple layers, the model is applied to each soil layer independently, starting at the upper layer. Vsai = Vsai-l + Rc -revap I: Model Development (T+ 17.8) (Tmx -~ } 0.6 mm (17) JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION 79 JAWRA Arnold, Srinivasan, Muttiah, and Williams where Tmx and Tmn are the daily maximum and minimum air temperatures in 'C. The model computes evaporation from soils and plants separately as described in Ritchie (1972), Potential soil water evaporation is estimated as a function of potential ET and leaf area index (area of plant leaves relative to the soil surface area). Actual soil water evaporation is estimated by using exponential functions of soil depth and water content. Plant water evaporation is simulated as a linear function of potential ET and leaf area index. Snow Melt- If snow is present, it may be melted on days when the second soil layer temperature exceeds O-C. Snow is melted as a function of the snow pack temperature using the equation SML = T (1.52 + 0.54 SPT) 0. ~ SML ~ SNO (18) where SML is the snowmelt rate in mm.d-l, SNO is the snow present in mm of water, T is the mean daily air temperature in °C,and SPT is the snow pack temperature in °C. The snow pack temperature is estimated with the equation SPT ==rnin {TB.T(2)} The weather variables for driving the hydrologic balance are precipitation, air tE)mperature, solar radiation, wind speed, and relative humidity. If daily precipitation and maximum/minimum temperature data are available, they can be input directly. If not, the weather generator can simulate daily rainfall and temperature. Solar radiation, wind speed, and relative humidity are always simulated. One set of weather variables may be simulated for the entire basin, or different weather may be simulated for each subbasin. Precipitation. The precipitation model developed by Nicks (1974) is a first-order Markov chain model. Thus, input to the model must include monthly probabilities of receiving precipitation if the previous day was dry and if the previous day was wet. Given the wet-dry state, the model determines stochastically if precipitation occurs or not. When a precipitation event occurs, the amount is determined by generating from a skewed normal daily precipitation distribution. The amount of daily precipitation is partitioned between rainfall and snowfall using average daily air temperature. (19) Air Temperature and Solar Radiation. Daily maximum and minimum air temperature and solar radiation are generated from a normal distribution corrected for wet-dry probability state. The correction factor is used to provide more deviation in temperatures and radiation when weather changes and for rainy days. Conversely, deviations are smaller on dry days. The correction factors are calculated to insure that long-term standard deviations of daily variables are maintained. where Ts is the temperature at the top of the snow pack and T(2) is the temperature at the center of soil layer 2. Melted snow is treated the same as rainfall for estimating runoff volume and percolation, but rainfall energy is set to 0.0 and peak runoff rate is estimated by assuming uniformly distributed rainfall for a 24-hour duration. Transmission Losses. Many semiarid watersheds have alluvial channels that abstract a considerable portion of streamflow (Lane, 1982). The abstractions, or transmission losses, reduce runoff volumes as the flood wave travels downstream. Lane's method described in USDA (1983) is used to estimate transmission losses. Channel losses are a function of channel width and length and flow duration. Both runoff and peak rate are adjusted when transmission losses occur. Wind Speed and Relative Humidity. Daily wind speed is simulated using a modified exponential equation given the mean monthly wind speed as input. The relative humidity model simulates daily average relative humidity from the monthly average by using a triangular distribution. As with temperature and radiation, the mean daily relative humidity is adjusted to account for wet- and dry-day effects. Ponds. Farm ponds are small structures that occur within a subbasin...Pond storage is simulated as a function of pond capacity, daily inflows and outflows, seepage, and evaporation. Ponds are assumed to have only emergency spillways. Required inputs are capacity and surface area. Surface area below capacity is estimated as a non-linear function of storage. JAWRA Weather Sedimentation Sediment Yield. Sediment yield is computed for each subbasin with the Modified Universal Soil Loss Equation (MUSLE) (Williams and Berndt, 1977). y= 80 11.8 JOURNAL (Vqp)O.56 (K) (C) OF THE AMERICAN (PE) WATER (LS) RESOURCES (20) ASSOCIATION ~ Large Area Hydrologic Modeling and Assessment -Part where y is the sediment yield from the subbasin in t, V is the surface runoff column for the subbasin in m3, qp is the peak flow rate for the subbasin in m3.s-1, K is the soil erodibility factor, C is the crop management factor, PE is the erosion control practice factor, and LS is the slope length and steepness factor. The LS factor is computed with the equation CWischmeier and Smith, 1978) Crop Growth Model The crop model is a simplification of the EPIC crop model (Williams et al., 1984). SWAT uses EPIC concepts of phenological crop development based on daily accumulated heat units, harvest index for partitioning grain yield, Monteith's approach (Monteith, 1977) for potential biomass, and water and temperature stress adjustments. A single model is used for simulating all the crops considered and SWAT is capable of simulating crop growth for both annual and perennial plants. Annual crops grow from planting date to harvest date or until the accumulated heat units equal the potential heat units for the crop. Perennial crops maintain their root systems throughout the year, although the plant may become dormant after frost. Phenological development of the crop is based on daily heat unit accumulation. It is computed using the equation (21) The exponent ~ varies with slope and is computed using the equation ~ = 0.6 [1 -exp(-35.835 8)] (22) The crop management factor, C, is evaluated for all days when runoff occurs using the equation, C = exp[(-0.2231 -CVM) I: Model Development exp(-O.OOl15 CV) + CVM HUi = ( Tmx.i +Tmn,i 2 ) -Tbj. HUh ?: 0 (25) (23) where HU, Tmx and Tmn are the values of heat units, maximum temperature, and minimum temperature in °C on day i and Tb is the crop-specific base temperature in °C(no growth occurs at or below Tb) of crop j. A heat unit index (HUI) ranging from O at planting to 1 at physiological maturity is computed as follows. where CM is the soil cover (above ground biomass + residue in kg.ha-l and CVM is the minimum value of C. The value of CVM is estimated from the average annual C factor using the equation (24) CVM = 1.4631n (CVA) + 0.1034 i The value of CVA for each crop is determined from tables prepared by Wischmeier and Smith (1978). Values of K are contained in the SCS Soils-5 database, and FE factors can be estimated for each subbasin using information contained in Wischmeier and Smith HUi Soil Temperature Daily average soil temperature is simulated at the center of each soil layer for use in hydrology and residue decay. The temperature of the soil surface is estimated using daily maximum and minimum air temperature and snow, plant, and residue cover for the day of interest plus the four days immediately preceding. Soil temperature is simulated for each layer using a function of damping depth, surface temperature, and mean annual air temperature. Damp~ ing depth is dependent upon bulk density and soil water. OF THE AMERICAN WATER RESOURCES ASSOCIATION (26) where HUI is the heat unit index for day i and PHU is the potential heat units required for maturity of crop j. The value of PHU is calculated by the model from normal planting and harvest dates. (1978). JOURNAL L HUh = h=l PHU j Potential Growth Interception of solar radiation is estimated with Beer's law equation (Monsi and Saeki, 1953) (27) where PAR is photosynthetic active radiation in MJ.m-2, RA is solar radiation in Ly, LAI is the leaf area index, and subscript i is the day of the year. 81 JAWRA ;t r 1; " ~ ~:: : ~ '111111 c II : .!;-:0 '1J \j Arnold, Srinivasan, Muttiah, and Williams "~ i ! t' j I} .,r - !;f :r " jr I,,! i .. Using Monteith's approach (Monteith, 1977), potential increase in biomass for a day can be estimated with the equation et al. (1976) and modified by Williams and (1978) for application to individual runoff events used to estimate organic N loss. The - i:i l if " 11 .:1 i:!, !i j, '" !C tJ.Bp,i=(BE;') (PARJ') estimates the daily organic N runoff loss based concentration of organic N in the top soil layer, sediment yield, and the enrichment ratio. Also, use of N is estimated using a supply and approach. The simulated nitrogen cycle is in Figure 4 and is a simplification of the actual nitrogen cycle. (28) where dBp is the daily potential increase in total biomass in kg.ha-l and BE is the crop parameter for converting energy to biomass in kg.m-2.ha-l.MJ. Potential biomass is adjusted daily due to stresses caused by water, nutrients, and temperature. LA! is simulated as a function of heat units and biomass. LAI is estimated with the equations LAI,- (LAImx}(BAG} +exp(9.5-.0006BAG} J -BAG , HIUi Phosphorus. The SWAT approach to estimating soluble p loss in surface runoff is based on the concept of partitioning pesticides into the solution and sediment phases as described by Leonard and Wauchope {Knisel, 1980). Because p is mostly associated with the sediment phase, the soluble p runoff is predicted using labile p concentration in the top soil layer, runoff volume, and a partitioning factor. Sediment , transport of p is simulated with a loading function as described in organic N transport. Crop use of p is I also estimated with the supply and demand approach. ~ Dl (29) 2 LAIi =(16) (LAImx) (l-HIUi) -~-~-, HIUi '~ > DLj (30) [1', ; t'l :ii jj , t where LAlmx is the maximum LAI potential for the crop, BAG is above ground biomass in kg.ha-l, and DLAI is the fraction of the growing season when LAI starts declining (=.75). The fraction of total biomass partitioned to the root system normally decreases from 0.3 to 0.5 in the seedling to 0.05 to 0.20 at maturity (Jones, 1985). The model estimates the root fraction to range linearly from 0.4 at emergence to 0.2 at maturity. Thus, the daily root fraction is computed with the equation RWTi = (0.4 -0.2 GLEAMS (Knisel, 1980) technology for simulating pesticide transport by runoff, percolate, soil evaporation, and sediment is used in the model (Figure 5). Pesticides may be applied at any time and rate to plant foliage or below the soil surface at any depth. The plant leaf-area-index determines what fraction of foliar applied pesticide reaches the soil surface. Also, a fraction of the application rate (called application efficiency) is lost to the atmosphere. Each pesticide has a unique set of parameters including solubility, half life in soil and on foliage, wash off fraction, organic carbon adsorption coefficient, and cost. Pesticide on plant foliage and in the soil degrade exponentially according to the appropriate half lives. Pesticide transported by water and sediment is calculated for each runoff event and pesticide leaching is estimated for each soil layer when percolation occurs. --, HUli) (31) where RWT is the fraction of total biomass partitioned to the root system on day i. Thus, BAG is calculated from the equation BAG = (1- RWTi) (BTO1i) (32) i\r~ ~~ where il, ':~: j!1 !. ,c. Nutr£ents Nitrogen. C:: '1 BTOT is total biomass Amounts in kg.ha-l on day i. of NO3-N contained Agricultural in runoff, The agricultural management component provides submodels that simulate tillage systems, application of irrigation water, fertilizer, and pesticides and grazing systems. lateral flow, and percolation are estimated as the products of the volume of water and the average concentration. Leaching and lateral subsurface flow in lower layers are treated with the same approach used in the upper layer, except that surface runoff is not considered. A loading function developed by McElroy JAWRA Management Tillage. The tillage component was designed to incorporate surface residue and chemicals into the soil. The user inputs the day of the tillage operation 82 JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION ~ Large Area Hydrologic Modeling and Assessment -Part I: Model Development NITROGEN Mineral N Organic N u N T"" I ~ :-:e 0) ~ 4) ~ ,Z ~ t) " ~ '.c ;.= z C) .~ bO .I ~ ~ ~ ~ b Q. h HUMUS Nitrification Figure 4. Nitrogen Forms and Transformations and selects the tillage implement from the database (over 100 implements are listed. Each implement has an associate mixing efficiency (0-1) which partitions the amount of residue that is incorporated into the soil and the remainder that is left on the surface. RSD = RSDo (1 -EF) Simulated with a SWAT Subbasin. where FC is the root zone field capacity (mm), SW is the root zone water content before irrigation. in mm, EFI is the efficiency ratio, and AIR is the volume of irrigation water applied ( mm). Fertilization. Fertilizer applications can also be scheduled by the user or automatically applied by the model. The user-scheduled option requires the user to input the application date, total amount of N and P, fraction of organic and inorganic N and P, and the soil layer of application. The model adds the amount of fertilizer to the proper nutrient pool (organic and inorganic) and to the specified soil layer. The automatic fertilization option requires the user to input the plant nitrogen stress level (0-1) to trigger fertilization, the amount of NO3 in the soil profile after fertilization, the soil layer of the application, the maximum amount of NO3 that can be applied in one year, and the minimum time between fertilizer applications. When the plant N stress level reaches the specified trigger level, the model automatically applies fertilizer to the NO3 storage of the specified soil layer to bring the entire profile to the specified level. (33) where RSD is the residue on the surface before tillage, RSDo is the residue after tillage, and EF is the mixing efficiency. Once the residue is incorporated, it has no impact on the model. Also, no adjustments are made to soil bulk density due to tillage. Irrigation. Both dryland and irrigated agriculture can be simulated. Irrigation applications may be scheduled by the user or automatically applied by the model. The user-scheduled option requires the user to input application dates, amounts, and application efficiencies. Irrigation water is added to fill the upper layer to field capacity and then added to fill the successive lower layers to field capacity until all of the water is applied. If automa.tic irrigation is specified, the user must input the application efficiency and a plant water stress level to trigger irrigation. When the user-specified stress level is reached, water is applied according to the equation . nl ANO3 = FNMX -L WNO3i (35) + ANO3 (36) i=l AIR= JOURNAL FC-SW 1- EFI OF THE AMERICAN (34) WATER RESOURCES ASSOCIATION WNO3fl 83 = WNO3fl JAWRA ~ ~ Arnold, Srinivasan, Muttiah, and Williams Figure 5. Pesticide Transfonnations Simulated Within a SWAT Subbasin. !:; f i j , ': :" I,I 'K 1 '1,c ~c' , 0'" where ANO3 is the amount of NO3 applied, FNMX is the amount of NO3 in the soil after fertilization, nl is the number of soil layers, and fl is the soil layer of the application. Organic N is also added according the amount of NO3 applied and the fraction of organic N (input). ON = ON + ANO3 ( torn 1- low, medium, or high level of p management, and the model automatically restores the upper two soil layers to 10, 20, and 30 ppm of labile P, respectively. Organic p is added like organic N in the previous equation. Pesticide Applications. The user inputs the pesticide number, the date, and the amount of pesticide applied. The user must also input the application efficiency factor to account for losses to the atmosphere. The amount of pesticide reaching the ground (added to the upper soil layer) and the amount intercepted by plants is computed as a function ofLAI. (37) torn where forn is the fraction of organic N in the total N application (0-1). Automatic p fertilization also occurs when the N stress level is reached. The user inputs a JAWRA 84 JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION Large Area Hydrologic Modeling and Assessment -Part with a modified form of the QUAL2E model (Ramanarayanan et al., 1996). The components include algae as chlorophyll-a dissolved oxygen, carbonaceous oxygen demand, organic nitrogen, ammonium nitrogen, nitrite nitrogen, nitrate nitrogen, organic phosphorus and soluble phosphorus. Water temperature is estimated from air temperature based on a relationship developed by Stefan and Preud'homme (1993) through regression analysis of numerous river observations. The relationship appears consistent with most rivers with the exceptions of natural springs and Grazing. Livestock grazing is simulated as a daily harvest operation. Users specify a daily grazing rate in kg.ha-l and the date grazing begins and ends. (38) BAG =BAG -BEAT where BEAT is the daily amount of biomass removed by livestock in kg.ha-l. Any number of grazing periods may occur during a year, and the grazing schedule may vary from year to year within a rotation. anthropological activity. Instream pesticide transformations are simulated with a modified form of a toxic model developed by Chapra (1989). The toxic is partitioned into dissolved and particulate in both the water and sediment layers. The major processes include reactions, volatilization, settling, diffusion, resuspension, and burial. Channel Routing. THe channel routing module consists of flood routing, sediment and chemical routing components. A ,detailed description of the routing components is found in Arnold et al. (1995). Channel Flood Routing. The flood routing model uses a variable storage coefficient method developed by Williams (1969). Channel inputs include the reach length, channel slope, bankfull width and depth, channel side slope, flood plain slope, and Manning's n for channel and floodplain. Flow rate and average velocity are calculated using Manning's equation and travel time is computed by dividing channel length by velocity. Outflow from a channel is also adjusted for transmission losses, evaporation, diversions, and return flow. Res~rvoir Routing Similar to channel routing, the reservoir routing module has water balance, sediment and chemical routing components. Reservoir Water Balance and Routing. The water balance for reservoirs includes inflow, outflow, rainfall on the surface, evaporation, seepage from the reservoir bottom, and diversions and return flow. There are currently three methods to estimate outflow. The first method simply reads in measured outflow and allows the model to simulate the other components of the water balance. The second method is for small uncontrolled reservoirs, and outflow occurs at a specified release rate when volume exceeds the principle storage. Volume exceeding the emergency spillway is released within one day. For larger managed reservoirs, a monthly target volume approach is used. Channel Sediment Routing. The sediment routing model (Arnold et at., 1995) consists of two components operating simultaneously (deposition and degradation). The deposition component is based on fall velocity and the degraqation component is based on Bagnold's stream power concept (Williams, 1980). Deposition in the channel and floodplain from the subbasin to the basin outlet is based on sediment particle fall velocity. Fall velocity is calculated as a function of particle diameter squared using Stokes Law. The depth of fall through a routing reach is the product of fall velocity and reach travel time. The delivery ratio is estimated for each particle size as a linear function of fall velocity, travel time, and flow depth. Stream power is used to predict degradation in the routing reaches. Bagnold (1977) defined stream power as the product of water density, flow rate, and water surface slope. Williams (1980) modified Bagnold's equation to place more weight on high values of stream power -stream power raised to 1.5. Available stream power is used to reentrain loose and deposited material until all of the material is removed. Excess stream power causes bed degradation. Bed degradation is adjusted by the USLE soil erodibility and cover factors Qf the channel and floodplain. Channel Instream JOURNAL Nutrient and Pesticide nutrient transformations are OF THE AMERICAN WATER RESOURCES Reservoir Sediment Routing. Inflow sediment yield to ponds and reservoirs (P/R) is computed with MUSLE. The outflow from P/R is calculated as the product of outflow volume and sediment concentration. Outflow P/R concentration is estimated using a simple continuity equation based on volumes and concentrations of inflow, outflow, and pond storage. Initial pond concentration is input and between storm concentration decreases as a function of time and median particle size of inflow sediment. Reservoir Nutrients and Pesticides. A simple model for phosphorus mass balance was taken from Thormann and Mueller (1987). The model assumes: (1) completely mixed lake; (2) phosphorus limited; and Routing. simulated ASSOCIATION I: Model Development 85 JAWRA Arnold, Srinivasan, Muttiah, and Williams 1!111 (3) total phosphorus can be a measure of trophic status. The first assumption ignores lake stratification and intensification of phytoplankton in the epilimnon. The second assumption is generally valid when nonpoint sources dominate and the third assumption implies that a relationship exists between total phosphorus and biomass. The phosphorus mass balance equation includes the concentration in the lake, inflow, outflow, and an overall loss rate. The lake toxic (pesticide) balance model is based on Chapra (1989) and assumes well mixed conditions. The system is partitioned into a well mixed surface water layer underlain by a well mixed sediment layer. The toxic is partitioned into dissolved and particulate in both the water and sediment layers. The major processes simulated by the model are loading, outflow, reactions, volatilization, settling, diffusion, resuspension, and burial. AIR= where AIR is the column of irrigation water to be applied, fd is an input parameter to allow for deficit irrigation, FC is the root zone field capacity, SW is the root zone soil water content before irrigation, and EFI is the irrigation efficiency (0-1). 4. Remove the water from the departure channel or reservoir. This is done by simply subtracting the actual transfer amount, or AIR, from the volume of water in the reservoir or the daily flow in the reach. MODEL LIMITATIONS SWAT is a long term water and sediment yield model that operates on a daily time step. Daily precipitation is input to the model and an empirical (curve number) equation is applied to daily rainfall without accounting for intensity. There are several reasons why the curve number approach was chosen over an infiltration equation for large area simulation that include: WATER TRANSFER AND MANAGEMENT For large basins it may be necessary to simulate water transfer. The transfer algorithm allows water to be transferred from any reach or reservoir to any other reach or reservoir in the watershed. It will also allow water to be diverted and applied directly to irrigate a subwatershed. There are four main steps to the algorithm: 1. Breakpoint rainfall (less than one day increments) is not readily available and is difficult to process. Storm disaggregation models have been developed (Obeysekera et at., 1987); however, they are stochastic with respect to intensities and often require inputs that are not readily available. 2. Subbasins are often relatively large (several km2) when simulating large river basins. It is relatively easy to obtain a weighted curve number and realistically simulate runoff. However, it is more difficult to "lump" saturated conductivity (a critical soil property for infiltration equations) since it can vary spatially by o;-ders of magnitude over relatively short distances. 3. Soils data is often not available with sufficient spatial detail for large basins to justify using an infiltration equation. 4. It relates runoff to soil type, land use, and management practices. 5. It is computation ally difficult. 1. Compute the maximum amount of water that can be transferred. This is the volume of water in the reservoir or the daily flow in the channel reach. 2. Compute the amount that is actually transferred. There are currently three options for determining the actual transfer amount: (a) specify the fraction of flow or volume to divert (0-1); (b) specify the minimum flow or volume remaining in the channel or reservoir after the water has been transferred; and (c) specify a daily amount to be diverted. More complex rules could be incorporated such as multiple destinations from multiple sources. An expert system may be appropriate for complex systems regarding order and amount to multiple destinations based on land use, previous weather conditions, soil water contents, reservoir levels, legal flow limits, etc. 3. Transfer the water to the destination. If the destination is a reach or reservoir, the actual transfer amount is added to the current storage in the reach or reservoir. If the destination is a subwatershed, a threshold must be reached before water is transferred and irrigation begins. If soil water content or crop stress drops below the input threshold, irrigation occurs. The amount of water needed for irrigation (Arnold and Stockle, 1992) is JAWRA Cfd) CFC-SW) lO-EFl One of the major limitacions to large area hydrologic modeling is the spatial variability associated with precipitation. There are more than 8000 raingage locations in the U .S. with more than 30 years of daily precipitation data. There are, on the average, two or three gages per county, which leaves several kilometers between gages. This can cause considerable errors in runoff estimation if one gage is used to 86 JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION Large Area Hydrologic Modeling and Assessment -J'art represent an entire subwatershed or even if an attempt is made to "spatially weight" precipitation for a subwatershed. Also, the data files are difficult to manipulate and contain considerable days of missing records. Weather generators can be extremely useful when measured data is unavailable and management scenarios are being compared. Daily weather generator parameters are available for generating weather sequences at a point. However, spatially correlated generators required for large area hydrologic simulation have not been developed. The physical processes driving large area weather phenomenon are not fully understood and many technical obstacles need to be overcome before spatially correlated rainfall generation is possible. Another possibility is to utilize the WSR-88D radar technology (formerly called NEXRAD -Next Generation Weather Radar) to measure aerial precipitation rates needed to drive large area hydrologic models. ARS researchers at Durant, Oklahoma, are currently testing WSR-88D and are simulating runoff based on WSR-88D estimates of precipitation. SWAT does not simulate detailed event-based flood and sediment routing. It was developed to predict agricultural management impacts on long-term (hundreds of years) erosion and sedimentation rates. The model operates on a daily time step, although a shorter and more flexible time increment would be a major enhancement to the model. The sediment routing equations are relatively simplistic and assume that the channel dimensions are static throughout the simulation. This may be unrealistic since simulations may be made for 100 years or more. The addition of algorithms to simulate channel downcutting and side slope stability would allow channel dimension to be continuously updated. Another limitation is the simplistic way the channel bed is described. The erodibility factor should be replaced with more detailed models that account for cohesive, noncohesive, and armored channels. Reservoir routing was originally developed for small reservoirs and assumes well-mixed conditions. The reservoir outflow calculations are simplistic and do not account for controlled operation. To adequately simulate large reservoirs, these items need to be addressed. SUMMARY time step and allows a basin to be subdivided into grid cells or natural subwatersheds. Major components of the hydrologic balance and their interactions are simulated including surface runoff, lateral flow in the soil profile, groundwater flow, evapotranspiration, channel routing, and pond and reservoir storage. The primary considerations in model development were to stress (1) land management, (2) water quality loadings, (3) flexibility in basin discretization, and (4) continuous time simulation. An attempt was made to simulate the major hydrologic components and their interactions as simply and yet as realistically as possible. An attempt was also made to use inputs that are readily available over large areas so the model can be used routinely in planning and decision making. SWAT is currently being utilized in several large area projects. SWAT provides the modeling capabilities of the HUMUS (Hydrologic Unit Model of the United States) project (Srinivasan et at., 1993). The HUMUS project simulates the hydrologic budget (Arnold et at., 1996) and sediment movement for approximately 2,100 8-digit hydrologic unit areas that have been delineated by the USGS. Findings of the project will be used in the Resource conservation Act (RCA) Assessment conducted by the Natural Resources Conservation Commission, scheduled for completion in 1991. Scenarios include projected ag and municipal water use, tillage and cropping system trends. The model is being used by NOAA to estimate nonpoint source loadings into all U. S. coastal areas as part of the National Coastal Pollutant Discharge Inventory. EPA is also utilizing the model to determine the severity of sediment contamination by pesticides in the U. S. OF THE AMERICAN Abbott, M. B., J. C. Bathurst, J. A. Cunge, P. E. O'Connell, and J. Rasmussen, 1986a. An Introduction to the European Hydrological System-Systeme Hydrologique European 'SHE' 1: History and Philosophy of a Physically-Based, Distribu ted Modeling System. J. Hydrol. 87:45-59. Abbott, M. B., J. C. Bathurst, J. A. Cunge, P. E. O'Connell, and J. Rasmussen, 1986b. An Introduction to the European Hydrological System-Systeme Hydrologique Europeen, 'SHE' 2: Structure of Physically-Based, Distributed Modeling System. J. Hydrol. 87:61-77. Arnold, J. G., P. M. Allen, and G. Bernhardt, 1993. A Comprehensive Surface-Groundwater Flow Model. J. Hydrol. 142:47-69. Arnold, J. G., R. Srinivasan, R. S. Muttiah, and P. M. Allen, 1996. Continental Scale Simulation of the Hydrologic Balance. ASCE J. ofHydrologic Eng. (accepted for publication). Arnold, J. G. and C. 0. Stockle, 1992. Simulation of Supplemental Irrigation from On-Fann Ponds. J. Irrigation and Drainage Div., ASCE 117(3):408-424. AND CONCLUSIONS WATER RESOURCES ASSOCIATION """ , "! c , :' i \' , i:;;: LITERATURE CITED A conceptual, continuous time model caned SWAT (Soil and Water Assessment Tool) was developed to assist water resource managers in assessing water supplies and nonpoint source pollution on watersheds and large river basins. The model operates on a daily JOURNAL I: Model Development 87 JAWRA J I Arnold, Srinivasan, Muttiah, and Williams I ~ " r\ : Arnold, J. G., J. R. Williams, and D. R. Maidment, 1995. Continuous-Time Water and Sediment Routing Model for Large Basins. ASCE J. of Hydr. Engr. 121(2). Arnold, J. G., J. R. Williams, A. D. Nicks, and N. B. Sammons, 1990. SWRRB: A Basin Scale Simulation Model for Soil and Water Resources Management. Texas A&M University Press, College Station, Texas. Bagnold, R. A., 1977. Bedload Transport in Natural Rivers. Water Resources Res. 13(2):303-312. Barros, A. P. and D. P. Lettenmaier, 1993. Dynamic Modeling of the Spatial Distribution ofPrecipitation in Remote Mountain Areas. Monthly Weather Review (121):1195-1214. Beasley, D. B., L. F. Huggins, and E. J. Monke, 1980. ANSWERS: A Model for Watershed Planning. Transactions of the ASAE. 23(4):938-944. Beven, K. J. and M. J. Kirk by, 1979. A Physically-Based Variable Contributing Area Model of Basin Hydrology. Hydrol. Sci. Bull. 24(1):43-69. Beven, K. J ., A. Calver, and E. M. Morris, 1987. The Institute of Hydrology Distributed Model (HDM). Report 98, Institute of Hydrology, Wallingford. United Kingdom. Beven, K. J., M. J. Kirkly, N. Schofield, and A. F. Tagg, 1984. Testing a Physically-Based Flood Forecasting Model (TOPMODEL) for Three UK Catchments. J. Hydrol. 69:119-143. Binley, A., J. Elgy, and K. J. Beven, 1989. A Physically Based Model of Heterogeneous Hill slopes. 1. Runoff Production. Water Resources Res. 25(6):1219-1226. Burnash, R. J. C., R. L. Ferral, and R. A. McGuire, 1973. A General Streamflow Simulation System- Conceptual Modeling for Digital Computers. Report by the Joint Federal State River Forecasts Center, Sacramento, California. Chapra, S. C., 1989. Water Quality Modeling of Toxic Organics in Lakes. CADSWES Working Paper No.4, University of Colorado, Boulder, Colorado. Crawford, N. H. and R. S. Linsley, 1966. Digital Simulation in Hydrology: The Stanford Watershed Model IV. Technical Report No, 39, Department of Civil Engineering, Stanford University, Palo Alto, California. Hargreaves, G. H. and Z. A. Samani, 1985. Reference Crop Evapotranspiration from Temperature. Applied Engr. Agric. 1:96.99. Hydrologic Engineering Center, 1981. HEC-1, Flood Hydrograph Package -Users Manual. U.S. Army Corps of Engineers, Davis, California. Jain, S. K., B. Storm, J. C. Bathurst, J. C. Refsgaard, and R. D. Singh, 1992. Application of the SHE to Catchments in India. Part 2. Field Experiments and Simulation Studies with the SHE on the Kolar Subcatchment of the Narmada River. Journal ofHydrology 140(1992):25-47. Jayatilaka, C. J., R. W. Gillham, D. W. Blowes, and R. J. Nathan, 1996. A Deterministic-Empirical Model of the Effect of the Capillary Fringe on Near-Stream Area Runoff. 2. Testing and Application. Journal ofHydrology 184(1996):317-336. Johansen, N. B., J. C. Imhoff, J. L. Kit tIe, and A. S. Donigian, 1984. Hydrological Simulation Program-Fortran (HSPF): User's Manual. EPA-600/3-84-066, U.S. EPA, Athens, Georgia. Jones, C. A., 1985. C-4 Grasses and Cereals. John Wiley and Sons, Inc., New York, New York, 419 pp. Knisel, W. G., 1980. CREAMS, A Field Scale Model for Chemicals, Runoff, and Erosion~from Agricultural Management Systems. U.S. Dept. Agric. Conserv. Res. Report No.26. Lane, L. J., 1982. Distributed Model for Small Semi-Arid Watersheds. ASCE J. Hydr. Div. 108(HY10):1114.1131. Laurenson, E. M. and R. G. Mein, 1983. RORB-Version 3: Runoff Routing Program -User Manual (2nd Edition). Department of Civil Engineering, Monash University, Monash, Australia. JAWRA 88 Leonard, R. A., W. G. Knisel, and D. A. Still, 1987. GLEAMS: Groundwater Loading Effects on Agricultural Management Systerns. Trans. ASAE 30(5):1403-1428. McElroy, A. D., S. Y. Chiu, J. W. Nebgen, A. Aletic, and R. W. Bennett, 1976. Loading Functions for Assessment of Water Pollution from Nonpoint Sources. Environ. Protection Tech. Serv., EPA 600/2-76-151. Monsi, M. and T. Saeki, 1953. Uber den Lictfaktor in den Pflanzengesellschaften und sein Bedeutung fur die Stoffproduktion.JapanJ.Bot. 14:22-52. Monteith, J. L., 1965. Climate and the Efficiency of Crop Production in Britain. Phil. Trans. Res. Soc. London Ser. B. 281:277329. Monteith, J. L., 1977. Climate and the Efficiency of Crop Production in Britain. Phil. Trans. Res. Sac. London Ser. B. 281:277329. Moore, R. J. and R. T. Clarke, 1981. A Distribution Function Approach to Rainfall-Runoff Modeling. Water Resour. Res., 17(5):1367-1382. Nicks, A. D., 1974. Stochastic Generation of the Occurrence, Pattern, and Location of Maximum Amount of Daily Rainfall. In: Proc. Symp. Statistical Hydrology, Aug.-Sept. 1971, Tucson, Arizona. U.S. Dept. Agric., Misc. Publ. No.1275, pp. 154-171. Obeysekera, J. T. B., G. Q. Tabios III, and J. D. Salas, 1987. On Parameter Estimation of Temporal Rainfall Models. Water Resources Res. 23(10):1837 -1850. Priestley, C. H. B. and R. J. Taylor, 1972. On the Assessment of Surface Heat Flux and Evaporation Using Large-Scale Parameters. Mon. Weather Rev. 100:81-92. Ramanarayanan, T. S., R. Srinivasan, and J. G. Arnold, 1996. Modeling Wister Lake Watershed Using a GIS-Linked Basin Scale Hydrologic/Water Quality Model. In: Proc. Third IntI. Conf. on Integrating Geographic Information Systems and Environmental Modeling. NCGIA, Santa Fe, New Mexico. Ritchie, J. T., 1972. A Model for Predicting Evaporation from a Row Crop with Incomplete Cover. Water Resources Res. 8:1204-1213. Rockwood, D. M., E. D. Davis, and J. A. Anderson, 1972. User Manual for COSSARR Model. U .S. Army Engineering Division, North Pacific, Portland, Oregon. Singh, V. P., 1989. Hydrologic Systems. VoI. 2, Watershed Modeling. Prentice Hall, Inc., Englewood Cliffs, New Jersey. Singh, V. P. (Editor), 1995. Computer Models of Watershed Hydrology. Chapter 1, Watershed Modeling. Water Resources Publications, Highlands Ranch, Colorado, pp. 1-22. Sloan, P. G., I. D. Morre, G. B. Coltharp, and J. D. Eigel, 1983. Modeling Surface and Subsurface Stormflow on Steeply-Sloping Forested Watersheds. Water Resouces Inst. Report 142, University ofKentucky, Lexington, Kentucky. Srinivasan, R., J. G. Arnold, R. S. Muttiah, C. Walker, and P. T. Dyke, 1993. Hydrologic Unit Model of the United States (HUMUS). In: Adv. in Hydro-Science and -Engineering, Sam Yang (Editor). Washington, D.C., June 7-11, Volume I, Part A, pp.451-456. Stefan, H. G. and E. B. Preud'homme, 1993. Stream Temperature Estimation from Air Temperature. Water Resources Bulletin 30(3):453-462. Sugawara, M., E. Ozaki, I. Wantanabe, and Y. Katsuyama, 1976. Tank Model and Its Application to Bird Creek, Wollombi Brook, Bihin River, Sanaga River, and Nam Mune. Research Note 11, National Center for Disaster Prevention, Tokyo, Japan. Thormann, R. V. and J. A. Mueller, 1987. Principles of Surface Water Quality Modeling and Control. Harper & Row, Inc., New York, New York. Todini, E., 1996. The ARNO Rainfall-Runoff Model. Journal of Hydrology 175(1996):339-382. JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION i Large Area Hydrologic Modeling and Assessment -Part I: Model Development u.s. Department of Agriculture, Soil Conservation Service, 1972. National Engineering Handbook, Hydrology Section 4, Chapters 4-10. U.S. Department of Agriculture, Soil Conservation Service, 1983. National Engineering Handbook, Hydrology Section 4, Chapter 19. U .S. Department of Agriculture, Soil Conservation Service, 1986. Urban Hydrology for Small Watersheds. Tech. Release 55. Williams, J. R., 1969. Flood Routing with Variable Travel Time or Variable Storage Coefficients. Trans. ASAE 12(1):100-103. Williams, J. R., 1980. SPNM, a Model for Predicting Sediment, Phosphorus, and Nitrogen from Agricultural Basins. Water Resources Bull. 16(5):843-848. Williams, J. R. and H. D. Berndt, 1977. Sediment Yield Prediction Based on Watershed Hydrology. Trans. ASAE 20(6):1100-1104. Williams, J. R. and R. W. Hann, 1973. HYMO: Problem-Oriented Language for Hydrologic Modeling -User's Manual. USDA, ARS-S-9. Williams,J. R. and R. W. Hann, 1978. Optimal Operation of Large Agricultural Watersheds with Water Quality Constraints. Texas Water Resources Institute, Texas A&M University, Tech. Report No.96. Williams, J. R., C. A. Jones, and P. T. Dyke, 1984. A Modeling Approach to Determining the Relationship Between Erosion and Soil Productivity. Trans. ASAE 27:129-144. Wischmeier, W. H. and D. D. Smith, 1978. Predicting Rainfall Erosion Losses: A Guide to Conservation Planning. U.S. Dept. Agric., Agric. Handbook No.537 . Young, R. A., C. A. Onstad, D. D. Bosch, and W. R. Anderson, 1987. AGNPS, Agricultural Non-Point Source Pollutin Model. A Watershed Analysis Tool. USDA Conservation Research Report 35. Zhao, R. J., 1984. Watershed Hydrological Modelling. Resources and Electric Power Press. Bejing, R.O.C. JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION ~, ~! Water 89 JAWRA