- Email: [email protected]
A system of erosion—sediment yield models
SOIL TECHNOLOGY Soil Technology
11 (1997)
43-55
A system of erosion-sediment yield models J.R. Williams *,l, J.G. Arnold 2 USDA, Agricultural
Research
Seruice, 808 East Blackland
Road, Temple, lX 76502, USA
1. Introduction The EPIC model (Williams et al., 1984) was developed in the early 1980’s to assess the effect of erosion on productivity. Since the 1985 RCA application, the model has been expanded and refined to allow simulation of many processes important in agricultural management (Sharpley and Williams, 1990). EPIC is a continuous simulation model that can be used to determine the effect of management strategies on agricultural production and soil and water resources. The drainage area considered by EPIC is generally a field sized area, up to 100 ha, where weather, soils and management systems are assumed to be homogeneous. The major components in EPIC are weather simulation, hydrology, erosion-sedimentation, nutrient cycling, pesticide fate, plant growth, soil temperature, tillage, economics and plant environment control. EPIC can be used to compare management systems and their effects on erosion and sediment yield. The management components that can be changed are crop rotations, tillage operations, irrigation scheduling, drainage, furrow diking, liming, grazing, manure handling and nutrient and pesticide application rates and timing. Simulator for water resources in rural basins (SWRRB) was developed to predict the effect of management decisions on water and sediment yields with reasonable accuracy for ungaged, rural basins (Arnold et al., 1990; Williams et al., 1985). The model was developed by modifying the CREAMS daily rainfall model (Knisel, 1980) for application to large, complex, rural basins. The major changes involved were (a) the model was expanded to allow simultaneous computations on several subwatersheds and (b) components were added to simulate weather, return flow, pond and reservoir storage, crop growth, transmission losses and sediment movement through ponds, reservoirs, streams and valleys. SWRRB operates on a daily time step and is efficient enough to run for many years (100 or more). Since the model is continuous in time, it can determine the * Corresponding author. ’ Hydraulic Engineer. 2 Agricultural Engineer. 00933-3630/97/$17X@ PII SO933-3630(96)00114-6
Copyright
0 1997 Elsevier
Science B.V. All rights reserved.
44
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
43-55
impacts of management such as crop rotations, planting and harvest dates, and chemical application dates and amounts. Basins can be subdivided into subareas based on differences in land use, soils, topography and climate. Sediment and associated chemicals are routed to the basin outlet. SWRRB has been validated on basins up to 500 km2 (Arnold and Williams, 1987). Recently, the SWAT (soil and water assessment tool) (Arnold et al., 1993) model was developed to model streamflow and sediment transport from large basins (several thousand km’). SWAT is a continuous time routing model that operates with a daily time step. The model is constructed to accept daily simulated or measured outflows from sub basins. Although the SWRRB model currently resides in SWAT to provide daily water and sediment yield estimates for sub basins, other models, like EPIC can be interfaced easily. Thus, it is possible to simulate water and sediment yield from some sub basins with SWRREl, some with EPIC, and to input measured values for others. The SWAT routing structure is based on that of the hydrologic model (HYMO) (Williams and Hann, 1972). Streamflow is routed through reaches using the variable storage coefficient routing method (Williams, 1975a; Williams, 1975b) and sediment is routed using a modification of the stream power (Bagnold, 1977) based method developed by Williams (1978) and modified by Arnold et al. (1990). Geographic information systems (GIS) are used to automate input development and to display spatially varying outputs. Besides providing river basin simulation capabilities, SWAT also serves as a shell to link the other Temple models. 2. Model description The components of the field scale EPIC model can be placed into 10 major divisions for discussion: hydrology, weather, erosion, nutrients, pesticide fate, soil temperature, plant growth, tillage, plant environment control and economics. A detailed description of the EPIC components was given by Williams et al. (1990). The components of SWRRB can be placed into eight major divisions: hydrology, weather, sedimentation, soil temperature, crop growth, nutrients, pesticide fate and agricultural management. A detailed description of the SWRRB components was given by Arnold et al. (1990). Only the erosion, sediment yield, sediment routing and closely related model components are described here. The closely related model components are described briefly beginning with those common to both EPIC and SWRRB and proceeding to those unique to each model. The SWAT streamflow and sediment routing components based on previous work (Williams and Harm, 1972; Williams, 1978) Arnold et al., 1992 are also described. 2.1. Erosion/
sediment yield related components
2.1.1. Components common to EPIC and SWRRB 2.1.1.2. h-&ace runofS. Surface runoff from daily rainfall is predicted using a procedure similar to the CREAMS runoff model, option one (Knisel, 1980; Williams and Nicks, 1982). Like the CREAMS model, runoff volume is estimated with a modification of the
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
43-55
45
SCS curve number method (USDA Soil Conservation Service, 1972). The curve number varies non linearly from the 1 (dry) condition at wilting point to the 3 (wet) condition at field capacity and approaches 100 at saturation. The EPIC model also includes a provision for estimating runoff from frozen soil. Peak runoff rate predictions are based on a modification of the Rational Formula. The runoff coefficient is calculated as the ratio of runoff volume to rainfall. The rainfall intensity during the watershed time of concentration is estimated for each storm as a function of total rainfall using a stochastic technique. The watershed time of concentration is estimated using Manning’s Formula considering both overland and channel flow. 2.1.1.2. Evapotranspiration. The model offers four options for estimating potential evaporation: Hargreaves and Samani (1985); Penman (1948); Priestley and Taylor (1972) and Penman-Monteith (Monteith, 1965). The Penman and Penman-Monteith methods require solar radiation, air temperature, wind speed and relative humidity as input. If wind speed, relative humidity and solar radiation data are not available, the Hargreaves or Priestley-Taylor methods provide options that give realistic results in most cases. The model computes soil and plant evaporation separately as described by Ritchie (1972). 2.1.2. Components unique to EPIC 2.1.2.1. Wind erosion. The Manhattan-Kansas wind erosion equation (Woodruff and Siddoway, 19651, was modified by Cole et al. (1982) for use in the EPIC model. The original equation computes average annual wind erosion as a function of soil erodibility, a climatic factor, soil ridge roughness, field length along the prevailing wind direction and vegetative cover. The main modification of the model was converting from annual to daily predictions to interface with EPIC. Two of the variables, the soil erodibility factor for wind erosion and the climatic factor, remain constant for each day of a year. The other variables, however, are subject to change from day to day. The ridge roughness is a function of a ridge height and ridge interval. .Field length along the prevailing wind direction is calculated by considering the field dimensions and orientation and the wind direction. The vegetative cover equivalent factor is simulated daily as a function of standing live biomass, standing dead residue and flat crop residue. Daily wind energy is estimated as a nonlinear function of daily wind velocity. 2.1.2.2. Crop growth model. A single model is used in EPIC for simulating all the crops considered (corn, grain sorghum, wheat, barley, oats, sunflower, soybean, alfalfa, cotton, peanuts, potatoes, durham wheat, winter peas, faba beans, rapeseed, sugarcane, sorghum hay, range grass, rice, casava, lentils and pine trees). Of course, each crop has unique values for the model parameters. Energy interception is estimated as a function of solar radiation and the crop’s leaf area index. The potential increase in biomass for a day is estimated as the product of intercepted energy and a crop parameter for converting energy to biomass. The leaf area index is simulated with equations dependent upon heat units, the maximum leaf area index for the crop, a crop parameter that initiates leaf area index decline and five stress factors.
46
J.R. Williams,
J.G. Arnold/Soil
Technology
II (1997)
43-55
Crop yield is estimated using the harvest index concept. Harvest index increases as a non linear function of heat units from zero at planting to the optimal value at maturity. The harvest index may be reduced by water stress during critical crop stages (usually between 30 and 90% of maturity). The fraction of daily biomass growth partitioned to roots is estimated to range linearly between two fractions specified for each crop at emergence and at maturity. Root weight in a soil layer is simulated as a function of plant water use within that layer. Root depth increases as a linear function of heat units and potential root zone depth. The potential biomass is adjusted daily if one of the plant stress factors is less than 1.0 using the product of the minimum stress factor and the potential biomass. The water-stress factor is the ratio of actual to potential plant evaporation. The temperature stress factor is computed with a function dependent upon the daily average temperature, the optimal temperature and the base temperature for the crop. The N and P stress factors are based on the ratio of accumulated plant N and P to the optimal values. The aeration stress factor is estimated as a function of soil water relative to porosity in the root zone. Roots are allowed to compensate for water deficits in certain layers by using more water in layers with adequate supplies. Compensation is governed by the minimum root growth stress factor (soil texture and bulk density, temperature and aluminum toxicity). The soil texture-bulk density relationship was developed by Jones (1983). 2.1.2.3. Tillage. The EPIC tillage component was designed to mix nutrients and crop residue within the plow depth, simulate the change in bulk density and convert standing residue to flat residue. Other functions of the tillage component include simulating ridge height and surface roughness. Tillage operations convert standing residue to flat residue using an exponential function of tillage depth and mixing efficiency. When a tillage operation is performed, a fraction of the material (equal the mixing efficiency) is mixed uniformly within the plow depth. Also, bulk density is reduced as a function of mixing efficiency, bulk density before tillage and undisturbed bulk density. After tillage, bulk density returns to the undisturbed value at a rate dependent upon infiltration, tillage depth and soil texture. 2.1.3. Components unique to SWRRB 2.2.3.1. Transmission losses. Many semiarid watersheds have alluvial channels that abstract large volumes of streamflow (Lane, 1982). The abstractions, or transmission losses, reduce runoff volumes as the flood wave travels downstream. SWRRB uses Lanes method described in ch. 19 of the SCS Hydrology Handbook (USDA Soil Conservation Service, 1983) to estimate transmission losses. Channel losses are a function of channel width and length and flow duration. Both runoff volume and peak rate are adjusted when transmission losses occur. 2.1.3.2. Ponds and reservoirs. Farm 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
J.R. Williams, J.G. Arnold/Soil
Technology I1 (1997) 43-55
41
area below capacity is estimated as a non linear function of storage. Reservoirs are treated similarly except they have emergency and principal spillways. Thus, required inputs include volume and surface area at both spillway elevations and the principal spillway releaserate. 2.1.3.3. Crop growth model. A single model is used in SWRRB for simulating all crops. Energy interception is estimatedas a function of solar radiation and the crop’s leaf area index. The potential increase in biomass for a day is estimated as the product of intercepted energy and a crop parameter for converting energy to biomass.The leaf area index is simulatedwith equationsdependentupon heat units. Crop yield is estimatedusing the harvest index concept. Harvest index increasesas a non linear function of heat units from zero at planting to the optimal value at maturity. The harvest index may be reduced by water stressduring critical crop stages(usually between 30 and 90% of maturity). The fraction of daily biomass growth partitioned to roots is estimated to range linearly from 0.4 at emergence to 0.2 at maturity. Root weight in a soil layer is simulated as a function of plant water use within that layer. Root depth increasesas a linear function of heat units and potential root zone depth. The potential biomassis adjusted daily if one of the plant stressfactors is less than 1.0 using the product of the minimum stress factor and the potential biomass. The water-stressfactor is the ratio of actual to potential plant evaporation. The temperature stressfactor is computed with a function dependentupon the daily average temperature, the optimal temperature and the basetemperature for the crop. 2.2. Erosion /sediment yield components 2.2.1. Componentsunique to EPIC The EPIC component for water induced erosion simulateserosion causedby rainfall and runoff and by irrigation (sprinkler and furrow). To simulate rainfall/runoff erosion, EPIC contains three equations: the USLE (Wischmeier and Smith, 1978), the MUSLE (Williams, 1975a; Williams, 1975b) and the Onstad-Foster modification of the USLE (Onstad and Foster, 1975). Only one of the equations (user specified) interacts with other EPIC components. The three equations are identical except for their energy components.The USLE dependsstrictly upon rainfall as an indicator of erosive energy. The MUSLE usesonly runoff variables to simulate erosion and sedimentyield. Runoff variables increasethe prediction accuracy, eliminate the need for a delivery ratio (used in the USLE to estimate sediment yield) and enable the equation to give single storm estimatesof sedimentyields. The USLE gives only annual estimates.The Onstad Foster equation contains a combination of the USLE and MUSLE energy factors. Thus, the water erosion model usesan equation of the form X==X(K)(CE)(PE)(LS)(ROKF) Xz=EI
for USLE
X = 1.586( Q . q,)0’56DAo.1z
for MUSLE
X ==0.646EI + 0.45( Q . q,)“‘33
for Onstand- Foster,
(1)
48
J. R. Williams,
J.G. Arnold/Soil
Technology
1 I I1 997) 43-55
where Y is the sediment yield (t/ha), K is the soil erodibility factor, CE is the crop management factor, PE is the erosion control practice factor, LS is the slope length and steepness factor, ROKF is the coarse fragment factor, EI is the rainfall energy factor, DA is the drainage area (ha), Q is the runoff volume (mm> and 4, is the peak runoff rate (mm/h). The PE value is determined initially by considering the conservation practices to be applied. The value of LS is calculated with the equation (Wischmeier and Smith, 1978). LS =
$0.1 t(65.41S2 + 4.563 + 0.065), (2) i 1 where S is the land surface slope (m/m), h is the slope length (m) and 5 is a parameter dependent upon slope. The value of 5 varies with slope and is estimated with the equation ,$=0.3S/[S+exp(-1.47+61.095)]
+0.2.
(3) The crop management factor is evaluated for all days when runoff occurs by using the equation CE=exp[(lnOJ
-lnCE,,,j)exp(
-1.15CV)
+lnCE,,,j]
(4)
where CE,,, is the minimum value of the crop management factor for crop j and CV is the soil cover (above ground biomass plus residue) (t/ha). The soil erodibility factor, K, is evaluated for the top soil layer at the start of each year of simulation with the equation
K = (0.2 + 0.3 exp( -0.256 SAN (1 - SIL/lOO))) . . . 1.0 . . . 1.0 -
( cL~~sIL)o’3
0.25C C + exp(3.72 - 2.95C) 0.7SNl SNl + exp( - 5.51 + 22.9SNI)
’
(5)
where SAN, SIL, CLA and C are the sand, silt, clay and organic carbon contents of the soil (%) and SNl = SAN/lOO. Eq. (5) allows K to vary from about 0.1 to 0.5. The first term gives low K values for soils with high coarse-sand contents and high values for soils with little sand. The fine sand content is estimated as the product of sand and silt divided by 100. The expression for coarse sand in the first term is simply the difference between sand and the estimated fine sand. The second term reduces K for soils that have high clay to silt ratios. The third term reduces K for soils with high organic carbon contents. The fourth term reduces K further for soils with extremely high sand contents (SAN > 70%). The runoff model supplies estimates of Q and qp. To estimate the daily rainfall energy in the absence of time-distributed rainfall, it is assumed that the rainfall rate is exponentially distributed: r, = rp exp( f/k),
(6)
J.R. Williams,
J.G. Arnold/Soil
Technology
11 11997) 43-55
49
where r is the rainfall rate at time t (mm/h), r-r is the peak rainfall rate (mm/h) and k is the decay constant (h). Eq. (6) contains no assumption about the sequence of rainfall rates (time distribution). The USLE energy equation in metric units is RE=R
12.1 +8.910gX
AR
(7)
where RE is the rainfall energy for water erosion equations and R is a rainfall amount (mm) during a time interval t (h). The energy equation can be expressed analytically as RE = 12.+df+ 0
8.9Jwr
log rdt.
0
(8)
Substituting Eq. (6) into Eq. (8) and integrating gives the equation for estimating daily rainfall energy RE = R[ 12.1 + 8.9(log rP - 0.434)],
(9
where R is the daily rainfall amount (mm). The rainfall energy factor, EI, is obtained by multiplying Eq. (9) by the maximum 0.5 h rainfall intensity (r,,s) and converting to the proper units: EI = R[ 12.1 + 8.9(log r-r - 0.434)] ( r0,5)/10-O0.
(10)
To compute values for rP, Eq. (6) is integrated to give R= (r,)(k)
(11)
and R,=R[l
-exp(-t/k)]
(12)
The value of R0,5 can be estimated by using u,,~: (13)
Ro.5 = ao.,R
where uo,s is the maximum daily amount of rainfall during 0.5 h divided by the total rainfall for the day. To determine the value of r,,, Eqs. (13) and (11) are substituted into Eq. (12) to give rp.= -2Rln(l
-ao.s).
(14)
Since rainfall rates vary seasonally, uo,s is evaluated for each month using Weather Service information (U.S. Department of Commerce, 1979). The frequency with which the maximum 0.5 h rainfall amount occurs is estimated using the Hazen plotting position equation (Hazen, 1930) F=27,
(15)
where F is the frequency with which the largest of a total of T events occurs. The total number of events for each month is the product of the number of years of record and the average number of rainfall events for the month. To estimate the mean value of u~,~, it is necessary to estimate the mean value of R,,. The value of R,, can be computed easily
50
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
43-55
if the maximum 0.5 h rainfall amounts are assumed to be exponentially distributed. From the exponential distribution, the expression for the mean 0.5 h rainfall amount is R 0.5.j
R 0.5F.j =
~
-1nq-j
(16)
where Ro,,,j is the mean maximum 0.5 h rainfall amount, R,,,, j is the maximum 0.5 h rainfall amount for frequency F and subscript j refers to the month. The mean ao,s is computed with the equation '0.5 a O.S,j =2 Ej
j
’
(17)
where R is the mean amount of rainfall for each event (average monthly rainfall/average number of days of rainfall) and subscript j refers to the month. Daily values of ao,5 are generated from a triangular distribution. The base of the triangular distribution is established by examining upper and lower limits of ao.5. The lower limit determined by a uniform rainfall rate gives a ,,s equal to 0.5/24 or 0.0208. The upper limit of ao,5 is set by considering a large rainfall event. In a large event, it is highly unlikely that all the rainfall occurs in 0.5 h (a = 1). The upper limit of ao,5 can be estimated by substituting a high value for r,, (250 mm/h is generally near the upper limit of rainfall intensity) into Eq. (12). ao.5u = exp( - 125/~),
(18)
where ao,5u is the upper limit of ao,5. The peak of the triangular distribution is set to aasj for month j. The coarse fragment factor is estimated with the equation (Simanton et al., 1984) ROKF = exp( - O.O3ROK),
(19)
where ROK is the percent of coarse fragments in the surface soil layer. 2.2.1.1. Irrigation. Erosion caused by applying irrigation water in furrows is estimated with MUSLE (Williams, 1975a; Williams, 1975b): y= 1.586(Q~~,)“~““(DA)o~‘“(K)(CE)(PE)(LS),
(20)
where the volume of runoff is estimated as the product of the irrigation volume applied and the irrigation runoff ratio. The peak runoff rate is estimated for each furrow by using Manning’s equation and assuming that the flow depth is 0.75 of the ridge height and that the furrow is triangular. 2.2.2. Components common to SWRRB and SWAT The approach to estimating sub-area sediment yield is similar to that of EPIC. The only differences are: (1) SWRRB uses only MUSLE, (2) the soil profile is static, the soil erodibility factor is input and held constant, (3) coarse fragments are not considered and (4) erosion by irrigation is not considered.
J.R. Williams,
J.G. Arnold/Soil
Technology
II (1997)
43-55
51
2.3. Sediment routing 2.3.1. Reach The sediment routing model consists of two components operating simultaneously (deposition and degradation). Deposition in the stream channel is based on the fall velocity of the sediment particles (Arnold et al., 1990). With a temperature of 22°C and a sediment density of 1.2 t/m3, Stokes’ Law for fall velocity becomes V,=411(d2),
(21) where V, is the fall velocity in m/h and d is the sediment particle diameter. The depth ( y,) that sediment of particle size d will fall during time, TT, is Yf = (Vf)G-?
(22)
where yf is the fall depth in m and TT is the travel time through a reach in h. The sedimentdelivery ratio (DR) through the reach is estimated with the equations DR= 1 -OS(y,)/d, DR=OS(d,)/y,
if yfd,,
(23) (24)
where dq is the depth of flow. Finally, deposition is calculated with the equation: DEP=SED,(l
-DR).
(25) Stream power is used to predict degradation in the routing reaches.Williams (1980) used Bagnold’s definition of stream power to develop a method for determining degradation in channels. Bagnold (1977) defined streampower, SP, with the equation SE’= gqs, )
(26)
where g is the density of the water, q is the flow rate and S, is the water surface slope. By applying stream power to bed-load predictions (Bagnold, 1977) and estimating model parameters(Williams, 19801,the equation for sedimentre-entrained, DEG,, is DEG, = bs,Y1-5(dur)(w)(d,S,V,)1’5,
(27)
where b,, is a parameter dependenton maximum streampower for the reach and V, is the velocity in the channel, dur is the duration of flow in h and w is the channel width in m. The parameter b,, can be estimatedwith the equation
where S, is the slope of the channel and the subscript mx refers to the maximum flow expected in the reach for extreme events. The value of q is assumedto equal some maximum rainfall intensity (250 mm/h) and b,, becomes b,, = (69.44g DA Ss)-o’5, where DA is the drainage area into the reach in km2.
(29)
J.R. Williams,
52
J.G. Arnold/Soil
Technology
II (1997) 43-55
All of the stream power is used for re-entrainment of loose and deposited material until all of the material has been removed. When this occurs, degradation of the bed material, DEG,, begins and is calculated by DEG, = KC DEG, ,
(30)
where K and C are MUSLE (Williams and Bemdt, 1977) factors for the stream channel. Total degradation, DEG, is the sum of the re-entrainment and bed degradation components. This amount is also allowed to be redeposited before reaching the basin outlet. DEG = ( DEG, + DEG, ) ( i - DR) .
(31)
Finally, the amount of sediment reaching the basin outlet, SED,,,, is SEDou, = SED,, - DEP + DEG, where SED,
(32)
is the sediment entering the reach.
2.3.2. Reservoir The sediment balance equation for the reservoirs SR, = SR,- 1 + SR,, - SR,,,
is
- SRnn,,
(33)
where SR i is the total sediment in the reservoir, SR,- r is the total sediment in the reservoir on the previous day, SR, is the incoming sediment, SR,,, is the sediment transported in the sediment outflow and SR,,, is the amount of sediment deposited in the reservoir. Computation of deposition is similar to the fall velocity algorithm used in the reach routing. DR=
1 -0.5(Y,)/d,
DR = 0.5( d,)/Y,
if Yr d,,
where d, is the depth of the reservoir, equations DEP, = SR,(l DEP, = 0 where SR,
(34)
- DR)
(SR,)(V,)
is the equilibrium
(SR,)(V,)
(35) Deposition,
DEP,,
is calculated
> SR,,
I SR,,,
with
the (36) (37)
sediment concentration for the reservoir.
2.4. Components unique to SWAT 2.4.1. Reach rating structure A new routing structure is required to convert the model into a distributed parameter model. Channel routing in SWRRB was performed from the subwatershed outlets directly to the basin outlet for simplicity. A more conventional reach routing structure is required to simulate large river basins. The reach routing structure of SWAT, similar to the WEPP watershed version (Lane and Nearing, 1989), HYMO (Williams and Harm, 1973) and ROT0 (Arnold et al., 1990), routes and addes flow down through the watershed through reaches and reservoirs. A set of commands are used to control routing and adding the flows down through the watershed. A routing command language was
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
developed and a list of the commands and corresponding the following table. Command comment subbasin
route
routres
transfer
add
routsub
recall
save
Code
Hyd. Location
Input 1
Hyd. Storage Location Hyd. Storage Location Hyd. Storage Location Departure Node Number Hyd. Storage Location Hyd. Storage Location Hyd. Storage Location file name Hyd. Storage Location file name
Subbasin No.
43-55
53
input parameters are given in
Input 2
Input 3
Input 4
Reach No.
Inflow Hyd.
Fraction in Channel
Reservoir No.
Inflow Hyd.
Dest. Node Number
Rule Code
Irrigation Transfer Code Subbasin No. Weather Destination icode
Inflow Hyd. 1
Inflow Hyd. 2
-
Flow Amount Transfered -
outflow Subbasin
Inflow Subbasin
-
-
-
finish 2.4.2. Groundwater jlow The system simulated by SWAT consists of four control volumes that include the: (1) surface, (2) soil profile or root zone, (3) shallow aquifer and (4) deep aquifer. The percolate from the soil profile is assumed to recharge the shallow aquifer. Once the water percolates to the deep aquifer it is lost from the simulated system and cannot return. The water balance for the shallow aquifer is SAQ, = SAQ,- 1 + Rc - revap -
qti -
percgw ,
(38)
where SAQ is the shallow aquifer storage, Rc is the recharge, revap is the water flow from the shallow aquifer back to the soil profile, qrf is the groundwater flow, percgw is
54
J.R. Williams,
J.G. Arnold/Soil
Technology
I1 (1997)
43-55
the percolate to the deep aquifer and i is the day with all units in m3 day-‘. Groundwater flow on day i is estimated by (Arnold et al., 1993) 41
=
4i-1
e
-aAr + Rc( 1 .O - eCaA’) ,
where CY is the constant of proportionality or the reduction factor. The relationship for water table height is (Arnold et al., 1993) hi = hi- 1 eCaA’ + &(l.O
- emmAr).
Sangrey et al. (1984) used an equation to estimate the delay time for return flow in a precipitation/groundwater response model. They utilized an exponential decay weighting function proposed by Venetis (1969)
RC, = (1.0 - e(-l-o/“))Rci
+ &l.o/S)Rci-lr
where 6 is delay time or drainage time of the aquifer (Sangrey et al., 1984). This equation will affect only the timing of the return flow and not the total volume.
References Arnold, J.G. and Williams, J.R., 1987. Validation of SWRRB: Simulator for water resources in rural basins. J. Water Res. Plann. Manage. ACSW, 113(2): 243-256. Arnold, J.G., Engel, B.A. and Srinivasan, R., 1993. A continuous time, grid cell watershed model. In: Application of Advanced Information Technologies for Management of Natural Resources. Sponsored by ASAE, June 17-19, Spokane, WA. Arnold, J.G., Williams, J.R., Griggs, R.H. and Sammons, N.B., 1992. SWRRBWQ-A basin scale model for assessing management impacts on water quality. Model Documentation and User Manual-Draft. USDAARS, Temple, TX, 195 pp. Arnold, J.G., Williams, J.R., Nicks, A.D. and Sammons, N.B., 1990. SWRRB-A basin scale simulation model for soil and water resources management. Texas A&M University Press, College Station, TX, 255 pp. Bagnold, R.A, 1977. Bed-load transport by natural rivers. Water Resourc. Res., 13(2): 303-312. Cole, G.W., Lyles, L. and Hagen, L.G., 1982. A simulation model of daily wind erosion soil loss. ASAE Pap., 82: 2575. Hargreaves, G.H. and Samani, Z.A., 1985. Reference crop evapotranspiration from temperature. Appl. Eng. Agric., 1: 96-99. Hazen, A, 1930. Flood flows, a study of frequencies and magnitudes. John Wiley & Sons, New York. Jones, CA, 1983. Effect of soil texture on critical bulk densities for root growth. Soil Sci. Sot. Am. J., 47: 1208-1211. Xnisel, W.G. (Editor), 1980. CREAMS: A field scale model for chemicals, runoff and erosion from agricultural management systems. USDA, Conservation Research Rep. No. 26, 643 pp. Lane, L.J., 1982. Distributed model for small semi-arid watersheds. J. Hydraul. Eng. ASCE, 109 (HYlO): 1114-1131. Lane, L.J. and Nearing, M.A. (Editor), 1989. USDA-Water erosion prediction project: Hillslope profile version. NSERL Rep. No. 2., NSERL, W. Lafayette, IN. Monte&h, J.L., 1965. Evaporation and environment. Symp. Sot. Exp. Biol., 19: 205-234. Onstad, C.A. and Foster, G.R., 1975. Erosion modeling on a watershed. Trans. ASAE 18: 288-292. Penman, H.L., 1948. Natural Evaporation from Open, Bare Soil and Grass. Proc. R. Sot. (London) A193: 120-145. Priestley, C.H.B. and Taylor, R.J., 1972. On the assessment of surface heat flux and evaporation using large scale parameters. Mon. Weather Rev., 100: 81-92.
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
43-55
55
Ritchie, J.T., 1972. A model for predicting evaporation from a row crop with incomplete cover. Water Resourc. Res., 8: 1204-1213. Sangrey, D.A., Harrop-Williams, K.O. and Klaiber, J.A., 1984. Predicting ground-water response to precipitation. AXE J. Geotech. Eng., llO(7): 957-975. Sharpley, A.N. and Williams J.R. (Editors), 1990. EPIC-Erosion/Productivity Impact Calculator. 1. Documentation. USDA Tech. Bull, 1768: 235 pp. Simanton, J.R., Rawitz, E. and Shirley, E.D., 1984. Effects of rock fragments on erosion of semi-arid rangeland soils. In: Erosion and Productivity of Soils Containing Rock Fragments. Soil Sci. Sot. A. Madison, WI, ch. 7, pp. 65-72. U.S. Department of Agriculture, Soil Conservation Service, 1972. Hydrology. In: National Engineering Handbook. USA Government Printing Office, ch. 4-10. U.S. Department of Agriculture, Soil Conservation Service, 1983. Hydrology. Sect. 4. National Engineering Handbook. USA Government Printing Office, ch. 19. U.S. Department of Commerce, 1979. Maximum short duration rainfall. National Summary, Climatic Data. Venetis, C., 1969. A study of the recession of unconfined aquifers. Bull. Int. Assoc. Sci. Hydrol., 19(4): 119-125. Williams, J.R., 1975a. HYMO flood routing. J. Hydrol., 26: 17-27. Williams, J.R., 1975b. Sediment-yield prediction with universal equation using runoff energy factor. In Present and Prospective Technology for Predicting Sediment Yield and Sources, USDA, ARS-S-40, pp. 224-252. Williams, J.R., 1978. A sediment yield routing model. In: Proc. Speciality Conf. Verification of Mathematical and Physical Models in Hydraulic Engineering. ASCE, College Park, MD, August 9- 11, pp. 662-670. Williams, J.R., 1980. SPNM, a model for predicting sediment, phosphorus and nitrogen yields from agricultural basins. Water Res. Bull. AWRA, 16(5): 843-848. Williams, J.R. and Bemdt, H.D., 1977. Determining the universal soil loss equation’s length-slope factor for watersheds. In: G.R. Foster (Editor), Soil Erosion: Prediction and Control, Proc. Natl. Conf. on Soil Erosion, May 24-26, 1976, Purdue Univ., West Lafayette, IN, SCSA Special Publ. No. 21, pp. 217-225. Williams, J.R. and Hann, Jr., R.W., 1972. HYMO, a problem-oriented computer language for building hydrologic models. Water Resourc. Res., 8: 79-86. Williams, J.R. and Harm, R.W., 1973. HYMO: Problem oriented computer language for hydrologic modeling. USDA .4RS-S-9, 76 pp. Williams, J.R. and Nicks, A.D., 1982. CREAMS hydrology model-option one. In: V.P. Singh (Editor) Applied Modeling Catchment Hydrology. Proc. Int. Symp. Rainfall-Runoff Modeling, May 18-21, 1981, Mississippi State, MS, p. 69086. Williams, J.R., Jones, C.A. and Dyke, P.T., 1984. A modeling approach to determining the relationship between erosion and soil productivity. Trans. ASAE 27(l): 129-144. Williams, J.R., Jones, C.A. and Dyke, P.T., 1990. The EPIC model. In: A.N. Sharpley and J.R. Williams (Editors), EPIC-Erosion/Productivity Impact Calculator: 1 Model Documentation. USDA Tech. Bull. No. 1768, ch. 2, pp. 3-02, 235 pp. Williams, J.R., Nicks, A.D. and Arnold, J.G., 1985. SWRRB, a simulator for water resources in rural basins. ASCE Hydr. J., 11 l(6): 970-986. Wischmeier, W.H. and Smith, D.D., 1978. Predicting rainfall erosion losses. Agriculture Handbook 537, USDA, SEA, 58 pp. Woodruff, N.P. and Siddoway, F.H., 1965. A wind erosion equation. Soil Sci. Sot. Am. Proc., 29: 602-608.
11 (1997)
43-55
A system of erosion-sediment yield models J.R. Williams *,l, J.G. Arnold 2 USDA, Agricultural
Research
Seruice, 808 East Blackland
Road, Temple, lX 76502, USA
1. Introduction The EPIC model (Williams et al., 1984) was developed in the early 1980’s to assess the effect of erosion on productivity. Since the 1985 RCA application, the model has been expanded and refined to allow simulation of many processes important in agricultural management (Sharpley and Williams, 1990). EPIC is a continuous simulation model that can be used to determine the effect of management strategies on agricultural production and soil and water resources. The drainage area considered by EPIC is generally a field sized area, up to 100 ha, where weather, soils and management systems are assumed to be homogeneous. The major components in EPIC are weather simulation, hydrology, erosion-sedimentation, nutrient cycling, pesticide fate, plant growth, soil temperature, tillage, economics and plant environment control. EPIC can be used to compare management systems and their effects on erosion and sediment yield. The management components that can be changed are crop rotations, tillage operations, irrigation scheduling, drainage, furrow diking, liming, grazing, manure handling and nutrient and pesticide application rates and timing. Simulator for water resources in rural basins (SWRRB) was developed to predict the effect of management decisions on water and sediment yields with reasonable accuracy for ungaged, rural basins (Arnold et al., 1990; Williams et al., 1985). The model was developed by modifying the CREAMS daily rainfall model (Knisel, 1980) for application to large, complex, rural basins. The major changes involved were (a) the model was expanded to allow simultaneous computations on several subwatersheds and (b) components were added to simulate weather, return flow, pond and reservoir storage, crop growth, transmission losses and sediment movement through ponds, reservoirs, streams and valleys. SWRRB operates on a daily time step and is efficient enough to run for many years (100 or more). Since the model is continuous in time, it can determine the * Corresponding author. ’ Hydraulic Engineer. 2 Agricultural Engineer. 00933-3630/97/$17X@ PII SO933-3630(96)00114-6
Copyright
0 1997 Elsevier
Science B.V. All rights reserved.
44
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
43-55
impacts of management such as crop rotations, planting and harvest dates, and chemical application dates and amounts. Basins can be subdivided into subareas based on differences in land use, soils, topography and climate. Sediment and associated chemicals are routed to the basin outlet. SWRRB has been validated on basins up to 500 km2 (Arnold and Williams, 1987). Recently, the SWAT (soil and water assessment tool) (Arnold et al., 1993) model was developed to model streamflow and sediment transport from large basins (several thousand km’). SWAT is a continuous time routing model that operates with a daily time step. The model is constructed to accept daily simulated or measured outflows from sub basins. Although the SWRRB model currently resides in SWAT to provide daily water and sediment yield estimates for sub basins, other models, like EPIC can be interfaced easily. Thus, it is possible to simulate water and sediment yield from some sub basins with SWRREl, some with EPIC, and to input measured values for others. The SWAT routing structure is based on that of the hydrologic model (HYMO) (Williams and Hann, 1972). Streamflow is routed through reaches using the variable storage coefficient routing method (Williams, 1975a; Williams, 1975b) and sediment is routed using a modification of the stream power (Bagnold, 1977) based method developed by Williams (1978) and modified by Arnold et al. (1990). Geographic information systems (GIS) are used to automate input development and to display spatially varying outputs. Besides providing river basin simulation capabilities, SWAT also serves as a shell to link the other Temple models. 2. Model description The components of the field scale EPIC model can be placed into 10 major divisions for discussion: hydrology, weather, erosion, nutrients, pesticide fate, soil temperature, plant growth, tillage, plant environment control and economics. A detailed description of the EPIC components was given by Williams et al. (1990). The components of SWRRB can be placed into eight major divisions: hydrology, weather, sedimentation, soil temperature, crop growth, nutrients, pesticide fate and agricultural management. A detailed description of the SWRRB components was given by Arnold et al. (1990). Only the erosion, sediment yield, sediment routing and closely related model components are described here. The closely related model components are described briefly beginning with those common to both EPIC and SWRRB and proceeding to those unique to each model. The SWAT streamflow and sediment routing components based on previous work (Williams and Harm, 1972; Williams, 1978) Arnold et al., 1992 are also described. 2.1. Erosion/
sediment yield related components
2.1.1. Components common to EPIC and SWRRB 2.1.1.2. h-&ace runofS. Surface runoff from daily rainfall is predicted using a procedure similar to the CREAMS runoff model, option one (Knisel, 1980; Williams and Nicks, 1982). Like the CREAMS model, runoff volume is estimated with a modification of the
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
43-55
45
SCS curve number method (USDA Soil Conservation Service, 1972). The curve number varies non linearly from the 1 (dry) condition at wilting point to the 3 (wet) condition at field capacity and approaches 100 at saturation. The EPIC model also includes a provision for estimating runoff from frozen soil. Peak runoff rate predictions are based on a modification of the Rational Formula. The runoff coefficient is calculated as the ratio of runoff volume to rainfall. The rainfall intensity during the watershed time of concentration is estimated for each storm as a function of total rainfall using a stochastic technique. The watershed time of concentration is estimated using Manning’s Formula considering both overland and channel flow. 2.1.1.2. Evapotranspiration. The model offers four options for estimating potential evaporation: Hargreaves and Samani (1985); Penman (1948); Priestley and Taylor (1972) and Penman-Monteith (Monteith, 1965). The Penman and Penman-Monteith methods require solar radiation, air temperature, wind speed and relative humidity as input. If wind speed, relative humidity and solar radiation data are not available, the Hargreaves or Priestley-Taylor methods provide options that give realistic results in most cases. The model computes soil and plant evaporation separately as described by Ritchie (1972). 2.1.2. Components unique to EPIC 2.1.2.1. Wind erosion. The Manhattan-Kansas wind erosion equation (Woodruff and Siddoway, 19651, was modified by Cole et al. (1982) for use in the EPIC model. The original equation computes average annual wind erosion as a function of soil erodibility, a climatic factor, soil ridge roughness, field length along the prevailing wind direction and vegetative cover. The main modification of the model was converting from annual to daily predictions to interface with EPIC. Two of the variables, the soil erodibility factor for wind erosion and the climatic factor, remain constant for each day of a year. The other variables, however, are subject to change from day to day. The ridge roughness is a function of a ridge height and ridge interval. .Field length along the prevailing wind direction is calculated by considering the field dimensions and orientation and the wind direction. The vegetative cover equivalent factor is simulated daily as a function of standing live biomass, standing dead residue and flat crop residue. Daily wind energy is estimated as a nonlinear function of daily wind velocity. 2.1.2.2. Crop growth model. A single model is used in EPIC for simulating all the crops considered (corn, grain sorghum, wheat, barley, oats, sunflower, soybean, alfalfa, cotton, peanuts, potatoes, durham wheat, winter peas, faba beans, rapeseed, sugarcane, sorghum hay, range grass, rice, casava, lentils and pine trees). Of course, each crop has unique values for the model parameters. Energy interception is estimated as a function of solar radiation and the crop’s leaf area index. The potential increase in biomass for a day is estimated as the product of intercepted energy and a crop parameter for converting energy to biomass. The leaf area index is simulated with equations dependent upon heat units, the maximum leaf area index for the crop, a crop parameter that initiates leaf area index decline and five stress factors.
46
J.R. Williams,
J.G. Arnold/Soil
Technology
II (1997)
43-55
Crop yield is estimated using the harvest index concept. Harvest index increases as a non linear function of heat units from zero at planting to the optimal value at maturity. The harvest index may be reduced by water stress during critical crop stages (usually between 30 and 90% of maturity). The fraction of daily biomass growth partitioned to roots is estimated to range linearly between two fractions specified for each crop at emergence and at maturity. Root weight in a soil layer is simulated as a function of plant water use within that layer. Root depth increases as a linear function of heat units and potential root zone depth. The potential biomass is adjusted daily if one of the plant stress factors is less than 1.0 using the product of the minimum stress factor and the potential biomass. The water-stress factor is the ratio of actual to potential plant evaporation. The temperature stress factor is computed with a function dependent upon the daily average temperature, the optimal temperature and the base temperature for the crop. The N and P stress factors are based on the ratio of accumulated plant N and P to the optimal values. The aeration stress factor is estimated as a function of soil water relative to porosity in the root zone. Roots are allowed to compensate for water deficits in certain layers by using more water in layers with adequate supplies. Compensation is governed by the minimum root growth stress factor (soil texture and bulk density, temperature and aluminum toxicity). The soil texture-bulk density relationship was developed by Jones (1983). 2.1.2.3. Tillage. The EPIC tillage component was designed to mix nutrients and crop residue within the plow depth, simulate the change in bulk density and convert standing residue to flat residue. Other functions of the tillage component include simulating ridge height and surface roughness. Tillage operations convert standing residue to flat residue using an exponential function of tillage depth and mixing efficiency. When a tillage operation is performed, a fraction of the material (equal the mixing efficiency) is mixed uniformly within the plow depth. Also, bulk density is reduced as a function of mixing efficiency, bulk density before tillage and undisturbed bulk density. After tillage, bulk density returns to the undisturbed value at a rate dependent upon infiltration, tillage depth and soil texture. 2.1.3. Components unique to SWRRB 2.2.3.1. Transmission losses. Many semiarid watersheds have alluvial channels that abstract large volumes of streamflow (Lane, 1982). The abstractions, or transmission losses, reduce runoff volumes as the flood wave travels downstream. SWRRB uses Lanes method described in ch. 19 of the SCS Hydrology Handbook (USDA Soil Conservation Service, 1983) to estimate transmission losses. Channel losses are a function of channel width and length and flow duration. Both runoff volume and peak rate are adjusted when transmission losses occur. 2.1.3.2. Ponds and reservoirs. Farm 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
J.R. Williams, J.G. Arnold/Soil
Technology I1 (1997) 43-55
41
area below capacity is estimated as a non linear function of storage. Reservoirs are treated similarly except they have emergency and principal spillways. Thus, required inputs include volume and surface area at both spillway elevations and the principal spillway releaserate. 2.1.3.3. Crop growth model. A single model is used in SWRRB for simulating all crops. Energy interception is estimatedas a function of solar radiation and the crop’s leaf area index. The potential increase in biomass for a day is estimated as the product of intercepted energy and a crop parameter for converting energy to biomass.The leaf area index is simulatedwith equationsdependentupon heat units. Crop yield is estimatedusing the harvest index concept. Harvest index increasesas a non linear function of heat units from zero at planting to the optimal value at maturity. The harvest index may be reduced by water stressduring critical crop stages(usually between 30 and 90% of maturity). The fraction of daily biomass growth partitioned to roots is estimated to range linearly from 0.4 at emergence to 0.2 at maturity. Root weight in a soil layer is simulated as a function of plant water use within that layer. Root depth increasesas a linear function of heat units and potential root zone depth. The potential biomassis adjusted daily if one of the plant stressfactors is less than 1.0 using the product of the minimum stress factor and the potential biomass. The water-stressfactor is the ratio of actual to potential plant evaporation. The temperature stressfactor is computed with a function dependentupon the daily average temperature, the optimal temperature and the basetemperature for the crop. 2.2. Erosion /sediment yield components 2.2.1. Componentsunique to EPIC The EPIC component for water induced erosion simulateserosion causedby rainfall and runoff and by irrigation (sprinkler and furrow). To simulate rainfall/runoff erosion, EPIC contains three equations: the USLE (Wischmeier and Smith, 1978), the MUSLE (Williams, 1975a; Williams, 1975b) and the Onstad-Foster modification of the USLE (Onstad and Foster, 1975). Only one of the equations (user specified) interacts with other EPIC components. The three equations are identical except for their energy components.The USLE dependsstrictly upon rainfall as an indicator of erosive energy. The MUSLE usesonly runoff variables to simulate erosion and sedimentyield. Runoff variables increasethe prediction accuracy, eliminate the need for a delivery ratio (used in the USLE to estimate sediment yield) and enable the equation to give single storm estimatesof sedimentyields. The USLE gives only annual estimates.The Onstad Foster equation contains a combination of the USLE and MUSLE energy factors. Thus, the water erosion model usesan equation of the form X==X(K)(CE)(PE)(LS)(ROKF) Xz=EI
for USLE
X = 1.586( Q . q,)0’56DAo.1z
for MUSLE
X ==0.646EI + 0.45( Q . q,)“‘33
for Onstand- Foster,
(1)
48
J. R. Williams,
J.G. Arnold/Soil
Technology
1 I I1 997) 43-55
where Y is the sediment yield (t/ha), K is the soil erodibility factor, CE is the crop management factor, PE is the erosion control practice factor, LS is the slope length and steepness factor, ROKF is the coarse fragment factor, EI is the rainfall energy factor, DA is the drainage area (ha), Q is the runoff volume (mm> and 4, is the peak runoff rate (mm/h). The PE value is determined initially by considering the conservation practices to be applied. The value of LS is calculated with the equation (Wischmeier and Smith, 1978). LS =
$0.1 t(65.41S2 + 4.563 + 0.065), (2) i 1 where S is the land surface slope (m/m), h is the slope length (m) and 5 is a parameter dependent upon slope. The value of 5 varies with slope and is estimated with the equation ,$=0.3S/[S+exp(-1.47+61.095)]
+0.2.
(3) The crop management factor is evaluated for all days when runoff occurs by using the equation CE=exp[(lnOJ
-lnCE,,,j)exp(
-1.15CV)
+lnCE,,,j]
(4)
where CE,,, is the minimum value of the crop management factor for crop j and CV is the soil cover (above ground biomass plus residue) (t/ha). The soil erodibility factor, K, is evaluated for the top soil layer at the start of each year of simulation with the equation
K = (0.2 + 0.3 exp( -0.256 SAN (1 - SIL/lOO))) . . . 1.0 . . . 1.0 -
( cL~~sIL)o’3
0.25C C + exp(3.72 - 2.95C) 0.7SNl SNl + exp( - 5.51 + 22.9SNI)
’
(5)
where SAN, SIL, CLA and C are the sand, silt, clay and organic carbon contents of the soil (%) and SNl = SAN/lOO. Eq. (5) allows K to vary from about 0.1 to 0.5. The first term gives low K values for soils with high coarse-sand contents and high values for soils with little sand. The fine sand content is estimated as the product of sand and silt divided by 100. The expression for coarse sand in the first term is simply the difference between sand and the estimated fine sand. The second term reduces K for soils that have high clay to silt ratios. The third term reduces K for soils with high organic carbon contents. The fourth term reduces K further for soils with extremely high sand contents (SAN > 70%). The runoff model supplies estimates of Q and qp. To estimate the daily rainfall energy in the absence of time-distributed rainfall, it is assumed that the rainfall rate is exponentially distributed: r, = rp exp( f/k),
(6)
J.R. Williams,
J.G. Arnold/Soil
Technology
11 11997) 43-55
49
where r is the rainfall rate at time t (mm/h), r-r is the peak rainfall rate (mm/h) and k is the decay constant (h). Eq. (6) contains no assumption about the sequence of rainfall rates (time distribution). The USLE energy equation in metric units is RE=R
12.1 +8.910gX
AR
(7)
where RE is the rainfall energy for water erosion equations and R is a rainfall amount (mm) during a time interval t (h). The energy equation can be expressed analytically as RE = 12.+df+ 0
8.9Jwr
log rdt.
0
(8)
Substituting Eq. (6) into Eq. (8) and integrating gives the equation for estimating daily rainfall energy RE = R[ 12.1 + 8.9(log rP - 0.434)],
(9
where R is the daily rainfall amount (mm). The rainfall energy factor, EI, is obtained by multiplying Eq. (9) by the maximum 0.5 h rainfall intensity (r,,s) and converting to the proper units: EI = R[ 12.1 + 8.9(log r-r - 0.434)] ( r0,5)/10-O0.
(10)
To compute values for rP, Eq. (6) is integrated to give R= (r,)(k)
(11)
and R,=R[l
-exp(-t/k)]
(12)
The value of R0,5 can be estimated by using u,,~: (13)
Ro.5 = ao.,R
where uo,s is the maximum daily amount of rainfall during 0.5 h divided by the total rainfall for the day. To determine the value of r,,, Eqs. (13) and (11) are substituted into Eq. (12) to give rp.= -2Rln(l
-ao.s).
(14)
Since rainfall rates vary seasonally, uo,s is evaluated for each month using Weather Service information (U.S. Department of Commerce, 1979). The frequency with which the maximum 0.5 h rainfall amount occurs is estimated using the Hazen plotting position equation (Hazen, 1930) F=27,
(15)
where F is the frequency with which the largest of a total of T events occurs. The total number of events for each month is the product of the number of years of record and the average number of rainfall events for the month. To estimate the mean value of u~,~, it is necessary to estimate the mean value of R,,. The value of R,, can be computed easily
50
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
43-55
if the maximum 0.5 h rainfall amounts are assumed to be exponentially distributed. From the exponential distribution, the expression for the mean 0.5 h rainfall amount is R 0.5.j
R 0.5F.j =
~
-1nq-j
(16)
where Ro,,,j is the mean maximum 0.5 h rainfall amount, R,,,, j is the maximum 0.5 h rainfall amount for frequency F and subscript j refers to the month. The mean ao,s is computed with the equation '0.5 a O.S,j =2 Ej
j
’
(17)
where R is the mean amount of rainfall for each event (average monthly rainfall/average number of days of rainfall) and subscript j refers to the month. Daily values of ao,5 are generated from a triangular distribution. The base of the triangular distribution is established by examining upper and lower limits of ao.5. The lower limit determined by a uniform rainfall rate gives a ,,s equal to 0.5/24 or 0.0208. The upper limit of ao,5 is set by considering a large rainfall event. In a large event, it is highly unlikely that all the rainfall occurs in 0.5 h (a = 1). The upper limit of ao,5 can be estimated by substituting a high value for r,, (250 mm/h is generally near the upper limit of rainfall intensity) into Eq. (12). ao.5u = exp( - 125/~),
(18)
where ao,5u is the upper limit of ao,5. The peak of the triangular distribution is set to aasj for month j. The coarse fragment factor is estimated with the equation (Simanton et al., 1984) ROKF = exp( - O.O3ROK),
(19)
where ROK is the percent of coarse fragments in the surface soil layer. 2.2.1.1. Irrigation. Erosion caused by applying irrigation water in furrows is estimated with MUSLE (Williams, 1975a; Williams, 1975b): y= 1.586(Q~~,)“~““(DA)o~‘“(K)(CE)(PE)(LS),
(20)
where the volume of runoff is estimated as the product of the irrigation volume applied and the irrigation runoff ratio. The peak runoff rate is estimated for each furrow by using Manning’s equation and assuming that the flow depth is 0.75 of the ridge height and that the furrow is triangular. 2.2.2. Components common to SWRRB and SWAT The approach to estimating sub-area sediment yield is similar to that of EPIC. The only differences are: (1) SWRRB uses only MUSLE, (2) the soil profile is static, the soil erodibility factor is input and held constant, (3) coarse fragments are not considered and (4) erosion by irrigation is not considered.
J.R. Williams,
J.G. Arnold/Soil
Technology
II (1997)
43-55
51
2.3. Sediment routing 2.3.1. Reach The sediment routing model consists of two components operating simultaneously (deposition and degradation). Deposition in the stream channel is based on the fall velocity of the sediment particles (Arnold et al., 1990). With a temperature of 22°C and a sediment density of 1.2 t/m3, Stokes’ Law for fall velocity becomes V,=411(d2),
(21) where V, is the fall velocity in m/h and d is the sediment particle diameter. The depth ( y,) that sediment of particle size d will fall during time, TT, is Yf = (Vf)G-?
(22)
where yf is the fall depth in m and TT is the travel time through a reach in h. The sedimentdelivery ratio (DR) through the reach is estimated with the equations DR= 1 -OS(y,)/d, DR=OS(d,)/y,
if yf
(23) (24)
where dq is the depth of flow. Finally, deposition is calculated with the equation: DEP=SED,(l
-DR).
(25) Stream power is used to predict degradation in the routing reaches.Williams (1980) used Bagnold’s definition of stream power to develop a method for determining degradation in channels. Bagnold (1977) defined streampower, SP, with the equation SE’= gqs, )
(26)
where g is the density of the water, q is the flow rate and S, is the water surface slope. By applying stream power to bed-load predictions (Bagnold, 1977) and estimating model parameters(Williams, 19801,the equation for sedimentre-entrained, DEG,, is DEG, = bs,Y1-5(dur)(w)(d,S,V,)1’5,
(27)
where b,, is a parameter dependenton maximum streampower for the reach and V, is the velocity in the channel, dur is the duration of flow in h and w is the channel width in m. The parameter b,, can be estimatedwith the equation
where S, is the slope of the channel and the subscript mx refers to the maximum flow expected in the reach for extreme events. The value of q is assumedto equal some maximum rainfall intensity (250 mm/h) and b,, becomes b,, = (69.44g DA Ss)-o’5, where DA is the drainage area into the reach in km2.
(29)
J.R. Williams,
52
J.G. Arnold/Soil
Technology
II (1997) 43-55
All of the stream power is used for re-entrainment of loose and deposited material until all of the material has been removed. When this occurs, degradation of the bed material, DEG,, begins and is calculated by DEG, = KC DEG, ,
(30)
where K and C are MUSLE (Williams and Bemdt, 1977) factors for the stream channel. Total degradation, DEG, is the sum of the re-entrainment and bed degradation components. This amount is also allowed to be redeposited before reaching the basin outlet. DEG = ( DEG, + DEG, ) ( i - DR) .
(31)
Finally, the amount of sediment reaching the basin outlet, SED,,,, is SEDou, = SED,, - DEP + DEG, where SED,
(32)
is the sediment entering the reach.
2.3.2. Reservoir The sediment balance equation for the reservoirs SR, = SR,- 1 + SR,, - SR,,,
is
- SRnn,,
(33)
where SR i is the total sediment in the reservoir, SR,- r is the total sediment in the reservoir on the previous day, SR, is the incoming sediment, SR,,, is the sediment transported in the sediment outflow and SR,,, is the amount of sediment deposited in the reservoir. Computation of deposition is similar to the fall velocity algorithm used in the reach routing. DR=
1 -0.5(Y,)/d,
DR = 0.5( d,)/Y,
if Yr
where d, is the depth of the reservoir, equations DEP, = SR,(l DEP, = 0 where SR,
(34)
- DR)
(SR,)(V,)
is the equilibrium
(SR,)(V,)
(35) Deposition,
DEP,,
is calculated
> SR,,
I SR,,,
with
the (36) (37)
sediment concentration for the reservoir.
2.4. Components unique to SWAT 2.4.1. Reach rating structure A new routing structure is required to convert the model into a distributed parameter model. Channel routing in SWRRB was performed from the subwatershed outlets directly to the basin outlet for simplicity. A more conventional reach routing structure is required to simulate large river basins. The reach routing structure of SWAT, similar to the WEPP watershed version (Lane and Nearing, 1989), HYMO (Williams and Harm, 1973) and ROT0 (Arnold et al., 1990), routes and addes flow down through the watershed through reaches and reservoirs. A set of commands are used to control routing and adding the flows down through the watershed. A routing command language was
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
developed and a list of the commands and corresponding the following table. Command comment subbasin
route
routres
transfer
add
routsub
recall
save
Code
Hyd. Location
Input 1
Hyd. Storage Location Hyd. Storage Location Hyd. Storage Location Departure Node Number Hyd. Storage Location Hyd. Storage Location Hyd. Storage Location file name Hyd. Storage Location file name
Subbasin No.
43-55
53
input parameters are given in
Input 2
Input 3
Input 4
Reach No.
Inflow Hyd.
Fraction in Channel
Reservoir No.
Inflow Hyd.
Dest. Node Number
Rule Code
Irrigation Transfer Code Subbasin No. Weather Destination icode
Inflow Hyd. 1
Inflow Hyd. 2
-
Flow Amount Transfered -
outflow Subbasin
Inflow Subbasin
-
-
-
finish 2.4.2. Groundwater jlow The system simulated by SWAT consists of four control volumes that include the: (1) surface, (2) soil profile or root zone, (3) shallow aquifer and (4) deep aquifer. The percolate from the soil profile is assumed to recharge the shallow aquifer. Once the water percolates to the deep aquifer it is lost from the simulated system and cannot return. The water balance for the shallow aquifer is SAQ, = SAQ,- 1 + Rc - revap -
qti -
percgw ,
(38)
where SAQ is the shallow aquifer storage, Rc is the recharge, revap is the water flow from the shallow aquifer back to the soil profile, qrf is the groundwater flow, percgw is
54
J.R. Williams,
J.G. Arnold/Soil
Technology
I1 (1997)
43-55
the percolate to the deep aquifer and i is the day with all units in m3 day-‘. Groundwater flow on day i is estimated by (Arnold et al., 1993) 41
=
4i-1
e
-aAr + Rc( 1 .O - eCaA’) ,
where CY is the constant of proportionality or the reduction factor. The relationship for water table height is (Arnold et al., 1993) hi = hi- 1 eCaA’ + &(l.O
- emmAr).
Sangrey et al. (1984) used an equation to estimate the delay time for return flow in a precipitation/groundwater response model. They utilized an exponential decay weighting function proposed by Venetis (1969)
RC, = (1.0 - e(-l-o/“))Rci
+ &l.o/S)Rci-lr
where 6 is delay time or drainage time of the aquifer (Sangrey et al., 1984). This equation will affect only the timing of the return flow and not the total volume.
References Arnold, J.G. and Williams, J.R., 1987. Validation of SWRRB: Simulator for water resources in rural basins. J. Water Res. Plann. Manage. ACSW, 113(2): 243-256. Arnold, J.G., Engel, B.A. and Srinivasan, R., 1993. A continuous time, grid cell watershed model. In: Application of Advanced Information Technologies for Management of Natural Resources. Sponsored by ASAE, June 17-19, Spokane, WA. Arnold, J.G., Williams, J.R., Griggs, R.H. and Sammons, N.B., 1992. SWRRBWQ-A basin scale model for assessing management impacts on water quality. Model Documentation and User Manual-Draft. USDAARS, Temple, TX, 195 pp. Arnold, J.G., Williams, J.R., Nicks, A.D. and Sammons, N.B., 1990. SWRRB-A basin scale simulation model for soil and water resources management. Texas A&M University Press, College Station, TX, 255 pp. Bagnold, R.A, 1977. Bed-load transport by natural rivers. Water Resourc. Res., 13(2): 303-312. Cole, G.W., Lyles, L. and Hagen, L.G., 1982. A simulation model of daily wind erosion soil loss. ASAE Pap., 82: 2575. Hargreaves, G.H. and Samani, Z.A., 1985. Reference crop evapotranspiration from temperature. Appl. Eng. Agric., 1: 96-99. Hazen, A, 1930. Flood flows, a study of frequencies and magnitudes. John Wiley & Sons, New York. Jones, CA, 1983. Effect of soil texture on critical bulk densities for root growth. Soil Sci. Sot. Am. J., 47: 1208-1211. Xnisel, W.G. (Editor), 1980. CREAMS: A field scale model for chemicals, runoff and erosion from agricultural management systems. USDA, Conservation Research Rep. No. 26, 643 pp. Lane, L.J., 1982. Distributed model for small semi-arid watersheds. J. Hydraul. Eng. ASCE, 109 (HYlO): 1114-1131. Lane, L.J. and Nearing, M.A. (Editor), 1989. USDA-Water erosion prediction project: Hillslope profile version. NSERL Rep. No. 2., NSERL, W. Lafayette, IN. Monte&h, J.L., 1965. Evaporation and environment. Symp. Sot. Exp. Biol., 19: 205-234. Onstad, C.A. and Foster, G.R., 1975. Erosion modeling on a watershed. Trans. ASAE 18: 288-292. Penman, H.L., 1948. Natural Evaporation from Open, Bare Soil and Grass. Proc. R. Sot. (London) A193: 120-145. Priestley, C.H.B. and Taylor, R.J., 1972. On the assessment of surface heat flux and evaporation using large scale parameters. Mon. Weather Rev., 100: 81-92.
J.R. Williams,
J.G. Arnold/Soil
Technology
11 (1997)
43-55
55
Ritchie, J.T., 1972. A model for predicting evaporation from a row crop with incomplete cover. Water Resourc. Res., 8: 1204-1213. Sangrey, D.A., Harrop-Williams, K.O. and Klaiber, J.A., 1984. Predicting ground-water response to precipitation. AXE J. Geotech. Eng., llO(7): 957-975. Sharpley, A.N. and Williams J.R. (Editors), 1990. EPIC-Erosion/Productivity Impact Calculator. 1. Documentation. USDA Tech. Bull, 1768: 235 pp. Simanton, J.R., Rawitz, E. and Shirley, E.D., 1984. Effects of rock fragments on erosion of semi-arid rangeland soils. In: Erosion and Productivity of Soils Containing Rock Fragments. Soil Sci. Sot. A. Madison, WI, ch. 7, pp. 65-72. U.S. Department of Agriculture, Soil Conservation Service, 1972. Hydrology. In: National Engineering Handbook. USA Government Printing Office, ch. 4-10. U.S. Department of Agriculture, Soil Conservation Service, 1983. Hydrology. Sect. 4. National Engineering Handbook. USA Government Printing Office, ch. 19. U.S. Department of Commerce, 1979. Maximum short duration rainfall. National Summary, Climatic Data. Venetis, C., 1969. A study of the recession of unconfined aquifers. Bull. Int. Assoc. Sci. Hydrol., 19(4): 119-125. Williams, J.R., 1975a. HYMO flood routing. J. Hydrol., 26: 17-27. Williams, J.R., 1975b. Sediment-yield prediction with universal equation using runoff energy factor. In Present and Prospective Technology for Predicting Sediment Yield and Sources, USDA, ARS-S-40, pp. 224-252. Williams, J.R., 1978. A sediment yield routing model. In: Proc. Speciality Conf. Verification of Mathematical and Physical Models in Hydraulic Engineering. ASCE, College Park, MD, August 9- 11, pp. 662-670. Williams, J.R., 1980. SPNM, a model for predicting sediment, phosphorus and nitrogen yields from agricultural basins. Water Res. Bull. AWRA, 16(5): 843-848. Williams, J.R. and Bemdt, H.D., 1977. Determining the universal soil loss equation’s length-slope factor for watersheds. In: G.R. Foster (Editor), Soil Erosion: Prediction and Control, Proc. Natl. Conf. on Soil Erosion, May 24-26, 1976, Purdue Univ., West Lafayette, IN, SCSA Special Publ. No. 21, pp. 217-225. Williams, J.R. and Hann, Jr., R.W., 1972. HYMO, a problem-oriented computer language for building hydrologic models. Water Resourc. Res., 8: 79-86. Williams, J.R. and Harm, R.W., 1973. HYMO: Problem oriented computer language for hydrologic modeling. USDA .4RS-S-9, 76 pp. Williams, J.R. and Nicks, A.D., 1982. CREAMS hydrology model-option one. In: V.P. Singh (Editor) Applied Modeling Catchment Hydrology. Proc. Int. Symp. Rainfall-Runoff Modeling, May 18-21, 1981, Mississippi State, MS, p. 69086. Williams, J.R., Jones, C.A. and Dyke, P.T., 1984. A modeling approach to determining the relationship between erosion and soil productivity. Trans. ASAE 27(l): 129-144. Williams, J.R., Jones, C.A. and Dyke, P.T., 1990. The EPIC model. In: A.N. Sharpley and J.R. Williams (Editors), EPIC-Erosion/Productivity Impact Calculator: 1 Model Documentation. USDA Tech. Bull. No. 1768, ch. 2, pp. 3-02, 235 pp. Williams, J.R., Nicks, A.D. and Arnold, J.G., 1985. SWRRB, a simulator for water resources in rural basins. ASCE Hydr. J., 11 l(6): 970-986. Wischmeier, W.H. and Smith, D.D., 1978. Predicting rainfall erosion losses. Agriculture Handbook 537, USDA, SEA, 58 pp. Woodruff, N.P. and Siddoway, F.H., 1965. A wind erosion equation. Soil Sci. Sot. Am. Proc., 29: 602-608.