Technical Bulletin Number 1768

Introduction 139 Vind erosion 139 General concepts 139 Modifications 141 Simulation results 146 Notations 149 References 151 The nutrient component of EPIC 152 Abstract 152 Introduction 152 Nutrient model compartments 152 Organic transformations 152 Inorganic transformations 154 Fertilizer application 155 Plant uptake 155 Data requirements 156 Model testing 161 Conclusions 163 Acknowledgments 163 References 164 111 8. Estimation of soil pH changes in EPIC Abstract Introduction Method Program Test results References Appendix 9. A sensitivity analysis of EPIC Abstract Introduction Analytical methods Results Discussion Guidance for EPIC users Notations Acknowledgments 10. Evaluation of EPIC using a dryland wheat-sorghum-fallow crop rotation 167 167 167 169 171 172 173 174 178 178 178 178 181 181 187 189 190 191 12. Evaluation of EPIC nutrient projections using soil profiles for virgin and cultivated lands 01 the same soil series 217 Abstract Introduction Soil management Evaluation of nutrient projections References 13. Demonstration and validation of crop grain yield simulation by EPIC Abstract Introduction Values for parameters Simulating corn grain yield response to soil depth Simulating grain yield of six crop species Simulating corn grain yield response to irrigation Overall conclusions References Abstract 191 Introduction 191 14. Perspectives Description of the validation test 192 Reference Model performance and discussion 196 Conclusions 203 References 204 217 217 217 218 219 220 220 220 222 227 228 232 232 234 235 235 11. Evaluation of EPIC using a sagebrush range site 206 Abstract 206 Introduction 206 Study area and methods 206 Results and discussion 209 Conclusions 215 References 216 IV 1.INTRODUCTION A.N. Sharpley and J.R. Villiams Accurate estimates of future soil productivity are essential in agricultural decision-making and planning from the field scale to the national level. Soil erosion reduces soil productivity, but the relationship between erosion and productivity has not been well defined. Until the relationship is adequately defined, selecting management strategies to maximize long-term crop production will be impossible. According to the Soil and Vater Resources Conservation Act (RCA)5 a report on the status of soil and water resources in the United States was required by 1985. One important aspect of these resources is the effect of erosion on long-term soil productivity. In 1981, the National Soil-Erosion/Soil- Productivity Research Planning Committee documented what was known about the problem, identified what additional knowledge was needed, and outlined a research approach for solving the problem (Villiams 1981). One of the most urgent needs outlined in the report was the development of a mathematical model for simulating erosion, crop production, and related processes. This model was envisioned to be used for determining the relationship between erosion and productivity for the United States. A national ARS erosion/productivity modeling team^ was organized and began developing the model in 1981, setting four goals in the development process. The model was to be (a) physically based and capable of simultaneously and realistically simulating the processes involved in erosion by using readily available inputs; (b) capable of simulating the processes as they would occur over hundreds of years, if necessary, because erosion can occur relatively slowly; (c) applicable to a wide range of soils, climates, and crops encountered in the United States; and (d) efficient, convenient to use, and capable of assessing the effects of management changes on erosion and soil productivity. The model developed, EPIC (Erosion-Productivity Impact Calculator), consists of (a) physically based components for simulating erosion, plant growth, and related processes and (b) economic components both for assessing the cost of erosion and for determining optimal management strategies. The model was developed in time to analyze the relationship between erosion and productivity for the RCA-mandated report. Beyond the analysis for the RCA report, EPIC should be a useful decision-making tool for determining optimal management ^J.R. Villiams, J.M. Shaffer, K.G. Renard, G.R. Foster, J.M. Laflen, L. Lyles, C.A. Onstad, A.N. Sharpley, A.D. Nicks, C.A. Jones, C.V. Richardson, P.T. Dyke, K.R. Cooley, and S.J. Smith. strategies from the farm to the national level. For example, EPIC is capable of dealing with decisions involving drainage, irrigation, water yield, erosion control (wind and water), weather, fertilizer and lime applications, pest control, planting dates, tillage, and crop residue management. As a research tool, EPIC can be used in developing, testing, and refining model components for various processes; sensitivity analysis to determine the importance of experimental variables and their interactions; and designing field experiments to obtain maximum information for the minimum cost. This publication details the EPIC model and its components and, in addition, provides information on several components and use of the model. REFERENCE National Soil-Erosion-Soil Productivity Research Planning Committee, USDA-ARS. 1981. Soil erosion effects on soil productivity: A research perspective. J. Soil Vater Conserv. 36:82-90. 2. THE EPIC MODEL J.R. Villiams, C.A. Jones, and P.T. Dyke ABSTRACT EPIC (Erosion/Productivity Impact Calculator) is a comprehensive model developed to determine the relationship between soil erosion and soil productivity throughout the united States. It continuously simulates the processes associated with erosion, using a daily time step and readily available inputs. Since erosion can occur relatively slowly, the model can simulate the process over hundreds of years if necessary. EPIC is generally applicable, computationally efficient, and capable of computing the effects of management changes on outputs. EPIC is composed of (a) physically based components for simulating erosion, plant growth, and related processes and (b) economic components for assessing the cost of erosion, determining optimal management strategies, etc. The EPIC physical components include hydrology, weather simulation, erosion-sedimentation, nutrient cycling, plant growth, tillage, and soil temperature. lODEL DESCRIPTION EPIC is a fairly comprehensive model, developed specifically for application to the erosion/productivity problem (National Soil Erosion-Soil Productivity Research Planning Committee USDA-ARS 1981; Villiams et al. 1985). Thus, computational efficiency and user convenience were important considerations in designing the model. The model can be run on a variety of mainframes and microcomputers. User convenience features are described in a separate publication. (Previous versions of EPIC have been described in limited detail by Villiams 1983; Villiams and Renard 1985; Villiams et al. 1983a, 1983b, 1984a, 1984b.) The drainage area considered by EPIC is generally small («1 ha) because soils and management are assumed to be spatially homogeneous. In the vertical direction, however, the model is capable of working with any variation in soil properties, the soil profile being divided into a maximum of 10 layers. Vhen erosion occurs, the second layer thickness is reduced by the amount of the eroded thickness, and the top layer properties are adjusted by interpolation (according to the distance the first layer is moved into the second layer). Vhen the second layer thickness becomes zero, the top layer is moved into the third layer, etc. The EPIC model consists of numerous component models that pertain to the following major aspects of the erosion/ productivity relationship: hydrology, weather, erosion, nutrients, plant growth, soil temperature, tillage, economics, and plant environment control. Descriptions of the EPIC components and the mathematic relationships used to simulate the processes involved follow. Hydrology Surface Runoff The runoff model simulates surface runoff volumes and peak runoff rates, given daily rainfall amounts. Runoff volume is estimated by using a modification of the Soil Conservation Service (SCS) curve number technique (U.S. Department of Agriculture, Soil Conservation Service 1972). The technique was selected for use because (a) it is reliable and has been used for many years in the united States; (b) it is computationally efficient; (c) the required inputs are generally available; and (d) it relates runoff to soil type, land use, and management practices. The use of readily available daily rainfall data is a particularly important attribute of the curve number technique because for many locations, rainfall data with time increments of less than 1 day are not available. Also, rainfall data manipulations and runoff computations are more efficient for data taken daily than at shorter intervals. Peak discharge rate is estimated by using a modification of the Rational formula. A stochastic element is introduced to the Rational equation to allow realistic simulation of peak discharge rates, given only daily rainfall and monthly rainfall intensity information. Runoff Volume Surface runoff is predicted for daily rainfall by using the SCS curve number equation (U.S. Department of Agriculture, Soil Conservation Service 1972) , ^ (il^)\ R > 0.2S p.^ R + 0.8s Q = 0.0, R < 0.2s where Q is the daily runoff, R is the daily rainfall, and s is a retention parameter (see "Notations" section). 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 (U.S. Department of Agriculture, Soil Conservation Service 1972) 100 s = 254 ( 1) [2.2] CN The constant, 254, in equation 2.2 gives s in millimeters. Thus, R and Q are also expressed in millimeters. CN2--the curve number for moisture condition 2, or average curve number--can be obtained easily for any area by using the SCS hydrology handbook (U.S. Department of Agriculture, Soil Conservation Service 1972). The handbook tables consider soils, land use, and management. Assuming that the handbook CN2 value is appropriate for a 57i slope, we developed the following equation for adjusting that value for other slopes. CN2S = I (CN3 - CN2) [1 - 2 exp(- 13.86 S)]+ CN2 [2.3] where CN2S is the handbook CN2 value adjusted for slope, CN3 is the curve number for moisture condition 3 (wet), and S is the average slope of the watershed. Values of CNi, the curve number for moisture condition 1 fdry), and CN3 corresponding to CN2 are also tabulated in the handbook. For computing purposes, CNi and CN3 were related to CN2 with the equations 20(100 - CN2) CNi = CN2 - — [2.4] 100 - CN2 + exp[2.533 - 0.0636(100 - CN2)] CN3 = CN2 exp[0.00673(100 - CN2)] [2.5] Fluctuations in soil water content cause the retention parameter to change according to the equation FFC s = si (1 ) [2.6] FFC + exp[wi - W2 (FFC)] where Si is the value of s associated with CNi, FFC is the fraction of field capacity, and wi and W2 are shape parameters. FCC is computed with the equation SV - VP FFC = [2.7] FC - VP where SV is the soil water content in the root zone, VP is the wilting point water content (1,500 kPa for many soils) and FC is the field capacity water content (33 kPa for many soils). Values for wi and W2 are obtained from a simultaneous solution of equation 2.6 according to the assumptions that s=S2 when FFC=0.5 and 8=83 when FFC-1.0: ' ^-^ -1.0] wi = In ^Q _ S3 + W2 [2.8] ^ s7 ^ W2 = 2.0 In Í ^'^ 1 1- ^2 - 0.5 - In 1.0 S3 s7 1.0 [2.9] where S3 is the CN3 retention parameter. Equations 2.8 and 2.9 assure that CNi corresponds with the wilting point and that the curve number cannot exceed 100. The FFC value obtained in equation 2.7 represents soil water uniformly distributed through the top 1.0 m of soil. Runoff estimates can be improved ii the depth distribution of soil water is known. For example, water distributed near the soil surface results in more runoff than the same volume of water uniformly distributed throughout the top meter of soil. Also, a soil surface associated with such a uniform distribution of soil water results in more runoff than a soil surface that is dry. Since EPIC estimates water content of each soil layer daily, the depth distribution is available. The effect of depth distribution on runoff is expressed in the depth weighting function H ¿=1 FFC^(- V- i) FFC* = I ¿=1 Z. < 1.0m [2.10] U-1 H where FFC* is the depth weighted FFC value for use in equation 2.6, Z is the depth (m) to the bottom of soil layer ¿, and M is the number of soil layers. Equation 2.10 performs two functions: (1) it reduces the influence of lower layers because FFC^ is divided by Z^ and (2) it gives proper weight to thick layers relative to thin layers because FFC is multiplied by the layer thickness. There is also a provision for estimating runoff from frozen soil. If the temperature in the second soil layer is less than OoC, the retention parameter is reduced by using the equation Sf = s [1 - exp(-0.00292 s) ] [2.11] where Sf is the retention parameter for frozen ground. Equation 2.11 increases runoff for frozen soils but allows significant infiltration when soils are dry. The final step in estimating the runoff volume is an attempt to account for uncertainty. The retention parameter or curve number estimate is based on land use, management, hydrologie soil group, land slope, and soil water content and distribution and is adjustable for frozen soil. However, many complex natural processes and artificial diversions that affect runoff are not accounted for in the model. Thus, the final curve number estimate is generated from a triangular distribution to account for this uncertain variation. The mean of the triangle is the best estimate of curve number based on using equations 2.10, 2.7, 2.6, 2.3, 2.2, and 2.11. The extremes are ±5 curve numbers from the mean. The generated curve number is substituted into equation 2.2 to estimate runoff with equation 2.1. Peak Runoff Rate Peak runoff rate predictions are based on a modification of the Rational formula qp = (P) (r) (A) / 360 [2.12] where q^ is the peak runoff rate (m^/s), /? is a runoff coefficient expressing the watershed infiltration characteristics, r is the rainfall intensity (mm/h) for the watershed's time of concentration, and A is the drainage area (ha). The runoff coefficient can be calculated for each storm if the amount of rainfall and runoff are known: p = [2.13] Since R is input and Q is computed with equation 2.1, p can be calculated directly. Rainfall intensity can be expressed with •4- V« /^ ■»• /■\l «x-4- ■» i^-r»f^r»i ■rx the relationship Rtc r [2.14] where Rtc is the amount of rainfall (mm) during the watershed's time of concentration, tc (h). The value of Rtc can be estimated by developing a relationship with total R. The Veather Service's TP-40 (Hershfield 1961) provides accumulated rainfall amounts for various durations and frequencies. Generally, Rtc and R24 (24-h duration is appropriate for the daily time step model) are proportional for various frequencies. Thus, Rtc = ÛR24 [2-15] where a is a dimensionless parameter that expresses the proportion of total rainfall that occurs during tc- The peak runoff equation is obtained by substituting equations 2.13, 2.14, and 2.15 into equation 2.12. (û)(q)(A) ^ ^ qp = —- [2.16] 360 (tc) The time of concentration can be estimated by adding the surface and channel flow times: i-CC + tes [2.17] where tec is the time of concentration for channel flow and tes is the time of concentration for surface flow (h). The tee can be computed by using the equation Le tee = — [2.18] where Le is the average channel flow length for the watershed (km) and Ve is the average channel velocity (m/s). The average channel flow length can be estimated by using the equation Le = y(L) (Lea) [2.19] where L is the channel length from the most distant point to the watershed outlet (km) and Lea is the distance along the channel to the watershed centroid (km). Average velocity can be estimated by using Manning's equation and assuming a trapezoidal channel with 2:1 side slopes and a 10:1 bottom width/depth ratio. Substitution of these estimated and assumed values gives L'cc - 7(1) (Lea) (n)Q-^^ 0.489 (qe)^'^^^)^'^^^ [2.20] where n is Manning's n, qe is the average flow rate (m^/s), and a is the average channel slope (m/m). Assuming that Lça=0.5L and that the average flow rate is about 6.35 mm/h and is a function of the square root of drainage area, yields the final equation for tee= (A) (<^) A similar approach is used to estimate tes"- tes = ^ [2-22] Vs where A is the surface slope length (m) and Vg is the surface flow velocity (m/s). Considering a strip 1 m wide down the sloping surface and applying Manning's equation gives Vs = o 2.23] where qs is the average surface flow rate and S is the land surface slope (m/m). Assuming that the average flow rate is about 6.35 mm/h and making substitutions into equations 2.22 and 2.23 to convert from m^/s to mm/h and from s to h, give the equation for estimating tcs^ 18 (S)"-*^ [2.24] Although some of the assumptions used in developing equations 2.21 and 2.24 may appear liberal, equation 2.17 generally gives satisfactory results for small homogeneous watersheds. Since equations 2.21 and 2.24 are based on hydraulic considerations, they are more reliable than purely empirical equations. To properly evaluate a, variation in rainfall patterns must be considered. For some short duration storms, most or all the rain occurs during tc causing a to approach its upper limit of 1.0. Other storms of uniform intensity cause a to approach a minimum value. All other patterns cause higher a values than the uniform pattern, because Rtc is greater than R24 for all patterns except the uniform. By substituting the products of intensity and time into equation 2.15, an expression for the minimum value of a, ûmn? is obtained: 24 Thus, a ranges within the limits tc — < 0 < 1.0 24 [2.25] Although confined between limits, the value of a is assigned with considerable uncertainty when only daily rainfall and simulated runoff amounts are given. Thus, a is generated from a gamma function with the base ranging from tc/24 to 1.0. The peak of the a distribution changes monthly because of seasonal differences in rainfall intensities. The Weather Service (U.S. Department of Commerce 1979) provides information on monthly maximum rainfall intensities that can be used to estimate the peak a for each month. Since the water erosion model estimates the maximum 0.5-h amount of each daily rainfall (a K)? these estimates are used in calculating a. Besides the convenience of avoiding double calculation, it is important to assure that a r and a are closely related for each storm. The relationship between a r and a can be obtained from TP-40 (Hershfield 1961) by fitting a log function to the 10-year frequency rainfall distribution: 6 where Rt is the rainfall amount (mm) for any time t. Re is the 6-h rainfall amount (mm), and b is a parameter used to fit the TP-40 relationship at any location. The value of a is computed with the equation Û = a . -^ [2.27] ^5 Details of the procedure for estimating a r are given in the water erosion section of this chapter. Percolation The EPIC percolation component uses a storage routing technique to simulate flow through soil layers. Flow from a soil layer occurs when soil water content exceeds field capacity. Vater drains from the layer until the storage returns to field capacity. The reduction in soil water is simulated with the routing equation SV^ = (SVQ^ - FC^) exp(-At / TT^) + FC^ [2.28] where SV and SVo are the soil water contents at the end and the start of time interval At (24 h) and TT is travel time through layer ¿ (h). 10 Thus, daily percolation can be computed by taking the difference between SV and SVo Q¿ = (SVQ^ - FC^) [1.0 - exp(-At / TT^) ] [2.29] where 0 is the percolation rate for layer ¿ (in mm/d). Travel time through a layer is computed with the linear storage equation PO. - FC. ^ ^ TT, = -I ^ [2.30] where PO is the porosity (mm], FC is field capacity (mm), and SC is saturated conductivity--that is, rate of water drainage through a saturated layer--(mm/h). The routing process is applied from the soil surface layer by layer through the deepest layer. Since the saturated conductivity of some layers may be much lower than that of others, the routing scheme can lead to an impossible situation (porosity of low saturated conductivity layers may be exceeded). For this reason, a back pass is executed from the bottom layer to the surface. If a layer's porosity is exceeded, the excess water is transferred to the layer above. This process continues through the top layer. Saturated conductivity may be input or estimated for each soil layer by using the equation 12.7 (100 - CLA.) (SS.) SC, = ^ ^ ^ — [2.31] ^ 100 - CLA^ + exp[11.45 - 0.097 (100 - CLA^) ] where CLA is the percentage of clay in soil layer ¿ and SS is the soil strength factor (described in the Growth Constraints section of this chapter). Percolation is also affected by freezing temperature. Vater can flow into a frozen layer but is not allowed to percolate from the layer. 11 Lateral Subsurface Flov Lateral subsurface flow is calculated simultaneously with percolation. The lateral flow function (similar to equation 2.29) is expressed in the equation {^ = (SVQ^ - FC^) [1.0 - exp(-1.0 / TTJP ] [2.32] where QR is the lateral flow rate for soil layer I (mm/d) and TTp^ is the lateral flow travel time (d). The lateral flow travel time is estimated for each soil layer by using the equation TT R^ 1000 (CLA^) (SS^) CLA^ + exp(10.047 - 0.148 CLA^) + 10 [2.33] Equations 2.29 and 2.32 must be solved simultaneously to avoid one process dominating the other, simply because the solution occurs first. Thus, an equation for the sum of percolation and lateral flow is written as h' 'I - (SVo^ - FC^) (l.O - exp(^) exp(^)) Taking the ratio of QR/0 and substituting the resulting equation 2.34 leads to the equation [2.34] '. into 0 + 0 1.0 - exp(^^)^ [1.0 - exp(:^)J (SVQ^ - FC^) (l.O - exp(^) exp(lM)| [2.35] Solving for 0 gives the final percolation equation (SVQ^ - FC^) (1.0 - exp(i|^) exp(^)j (l.O - exp(^)) 2.0 exp(lM) _ exp(lM) [2.36] The calculated 0 value is substituted into equation 2.34 to obtain the final estimate of QR. 12 Evapotranspiration The model offers two options for estimating potential evaporation--the Priestley-Taylor (1972) and Penman (1948) methods. The Penman method requires solar radiation, air temperature, wind speed, and relative humidity as inputs. If wind speed and relative humidity data are not available, the Priestley-Taylor method provides an option that usually gives realistic results. The model computes evaporation from soils and plants separately, as described by Ritchie (1972). Potential soil water evaporation is estimated as a function of potential evaporation and leaf area index (LAI, 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 evaporation and leaf area index. Potential Evaporation The Penman (1948) option for estimating potential evaporation is based on the equation [2.37] \ - '7—' <^' * *7—' *'" *'» ■ '1' " 5+7 HV 0 + 7 where EQ is the potential evaporation (mm), 6 is the slope of the saturation vapor pressure curve (kPa/oC), 7 is a psychrometer constant (kPa/oC), ho is the net radiation (MJ/m2), G is the soil heat flux (MJ/m2), HV is the latent heat of vaporization (MJ/kg), f(V) is a wind speed function (mm/d/kPa), ea is the saturation vapor pressure at mean air temperature (kPa), and ed is the vapor pressure at mean air temperature (kPa). The latent heat of vaporization is estimated with the temperature function HV = 2.50 - 0.0022 T [2-38] where T is the mean daily air temperature (oC). The saturation vapor pressure is also estimated as a function of temperature by using the equation ea = 0.1 exp Í54.88 - 5.03 ln(T + 273) - ^"^^^ ] [2.39] M T + 273^ The vapor pressure is simulated as a function of the saturation value and the relative humidity: ed = (ea) (RH) [2-40] 13 where RH is the relative humidity expressed as a fraction. The slope of the saturation vapor pressure curve is estimated with the equation 6 = ( ^ ( —Ö-M- . 5.03) [2.41] T + 273 T + 273 The psychrometer constant is computed with the equation 7 = 6.6 X 10-4 PB [2.42] where PB is the barometric pressure (kPa). The barometric pressure is estimated as a function of elevation by using the equation PB = 101 - 0.0115 ELEV + 5.44 X 10-7 ELEV2 [2.43] where ELEV is the elevation of the site (m). The soil heat flux is estimated by using air temperature on the day of interest plus 3 days prior. G = 0.12(T. - (-i^i ^ ^)j [2.44] where T is the mean daily air temperature on day i (oC). Solar radiation is adjusted to obtain net radiation by using the equation 0.9 RA. h . = RA. (1.0 - AB.) - RAB. ( i + 0.1) [2.45] ^ ^ RAMX.1 where RA is the solar radiation (MJ/m2), AB is albedo, RAB is the net outgoing long wave radiation (MJ/m2) for clear days, and RAMX is the maximum solar radiation possible (MJ/m2) for the location on day i. The RAB value is estimated with the equation RAB. = 4.9 X 10-9 (0.34 - 0.14 Ve^ ) (T. + 273)^ [2.46] The maximum possible solar radiation is computed with the equations RAMX = 30 1.0 + 0.0335 sin[ — (i + 88.2) 1 365 14 XT sin(— LAT) sin(SD) + cos(— LAT) cos(SDrsin(XT) 360 360 [2.47] XT = cos-i (-tan( — LAT) tan(SD)] , 0 < XT < T [2.48] ^ 360 ^ where LAT is the latitude of the site in degrees, SD is the sun's declination angle (radians), and i is the day of the year. The sun's declination angle is calculated with the equation SD. = 0.4102 sin[ ^ (i - 80.25) ] [2.49] ^ 365 Finally, the wind function of the Penman equation is approximated with the relationship f(Y) = 2.7 +1.63 V [2.50] where V is the mean daily wind speed at a 10-m height (m/s). The Priestley-Taylor (1972) method provides estimates of potential evaporation without wind and relative humidity inputs. The simplified equation based only on temperature and radiation is E = 30.6 (h ) ( ^ ) [2.51] ° ^ 6 + 0.68 The neo radiation is estimated with the equation "oi = ¡í¡ "*! (1 - "i) P.62] instead of equation 2.45, which is used in the Penman method. Similarly, equation 2.41 is replaced to estimate the slope of the saturation vapor pressure curve with the equation s -_ exp Í21.3 - ^304_ ][ 5304 | ^3 53^ ^ (T + 273) ^HT + 273)2 f 15 Both methods estimate albedo by considering the soil, crop, and snow cover. If a snow cover exists with 5 mm or greater water content, the value of albedo is set to 0.6. If the snow cover is less than 5 mm and no crop is growing, the soil albedo is the appropriate value. Vhen crops are growing, albedo is determined by using the equation AB = 0.23 (1.0 - EA) + (ABs) (EA) [2.54] where 0.23 is the albedo for plants, ABg is the soil albedo, and EA is a soil cover index. The value of EA ranges from 0 to 1.0 according to the equation EA = exp(-0.1 CV) [2.55] where CV is the sum of the above ground biomass and crop residue (t/ha). Soil and Plant Evaporation The model computes evaporation from soils and plants separately by an approach similar to that of Ritchie (1972). Potential plant water evaporation is computed with the equations 0 < LAI < 3.0 [2.56] LAI > 3.0 [2.57] where Ep is the predicted plant water evaporation rate (mm/d). If soil water is limited, plant water evaporation will be reduced as described in the plant growth section of this chapter. Potential soil water evaporation is simulated by considering soil cover according to the following equation Eg = min[ (E^) (EA), E^ - Ep ] [2.58] where Eg is the potential soil water evaporation rate (mm/d). Actual soil water evaporation is estimated on the basis of the top 0.2 m of soil and snow cover, if any. If snow is present, it is evaporated at the potential soil water evaporation rate. Vhen all snow is evaporated, soil water evaporation begins. Such evaporation is governed by soil depth and water content according to the equation F ^v .(LAI) fcp - 3 .0 if- ^0' ^h - \ Z / 0.2 Z / 0.2 + exp[-2.92 - 1.43 (—) ] 0.2 [2.59] 16 where EV is the total soil water evaporation (mm) from soil of depth Z (m). Potential soil water evaporation for a layer is estimated by taking the difference between EV's at the layer boundaries: SEV^ = ^\¿) - ^\¿.,^ [2.60] where SEV is the potential soil evaporation for layer ¿ (mm). The depth distributed estimate of soil water evaporation may be reduced according to the following equation if soil water is limited in a layer: * (2.5 (SV. - FC.)) SEV^ = SEV^ exp —\ , SV^ < FC^ [2.61] FC^ - VP¿ where SEV/, is the adjusted soil water evaporation estimate (mm). SEV¡ = SEV^ , SV^ > FC^ [2.62] The final step in adjusting the evaporation estimate is to assure that the soil water supply is adequate to meet the demand: SEvJ = min(SEV^ , SV^ - 0.5 VP^) [2.63] Equation 2.63 allows soil in the top 0.2 m to dry to half the soil water content corresponding to the wilting point. Snovmelt The EPIC snowmelt component is similar to that of the CREAMS model (Knisel 1980). If snow is present, it is melted on days when the maximum temperature exceeds O.O^C by using the equations SML = 4.57 T^, SML < SNO [2.64] SML = SNO [2.65] where SML is the snowmelt rate (mm/d), Tmx is the daily maximum air temperature (^C), and SNO is the water content of snow before melt occurs (mm). 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-h duration. 17 Vater Table Dynamics The water table height is simulated without direct linkage to other soil water processes in the root zone to allow for offsite water effects. The model drives the water table up and down between input values of maximum and minimum depths from the surface. The driving mechanism is a function of rainfall, surface runoff, and potential evaporation, as given in the equation VTBL. - VTBL.^ - Vl(VTBL._j - VTL) [2.66] where VTBL is the depth (m) from the surface to the water table on day i, VI is the driving function, and VTL is the appropriate limit. The driving equations are VI = mii](0.1, I V2 I ) [2.67] ^2 ^ RFS - qS - EOS [2.68] EOS where RFS, QS, and EOS are the sums of rainfall, runoff, and potential evaporation for 30 days before day i and V2 is a scaling factor. Equation 2.68 causes the water table to rise faster than it falls because the denominator is larger during recession. The maximum water table depth, VTMX, is substituted into equation 2.66 for VTL when the water table is falling. Conversely, VTL is set to the minimum water table depth, VTMN, on the rising side. VTL = VTMX , V2 < 0.0 [2.69] VTL = VTMN , V2 > 0.0 [2.70] Obviously, equation 2.66 gives highest rising rates when V2 is large and when VTBL«VTMX. As VTBLKVTMN the rate of rise approaches zero. The reverse is true on the falling side. y ., The weather variables necessary for driving the EPIC model are weat er precipitation, air temperature, and solar radiation. If the Penman method is used to estimate potential evaporation, wind speed and relative humidity are also required. Of course, wind speed is also needed when wind-induced erosion is simulated. If daily precipitation, air temperature, and solar radiation data are available, they can be input directly into EPIC. Rainfall and temperature data are available for many areas of the united States, but solar radiation, relative humidity, and wind data are scarce. Even rainfall and temperature data are generally not 18 adequate for the long-term EPIC simulations (100 years+). Thus, EPIC provides options for simulating various combinations of the five weather variables. Descriptions of the models used for simulating precipitation, temperature, radiation, relative humidity, and wind follow. Precipitation The EPIC precipitation model developed by Nicks (1974) is a first-order Markov chain model. Thus, input for the model must include monthly probabilities of receiving precipitation. On any given day, the input must include information as to whether the previous day was dry or wet. A random number (0-1) is generated and compared with the appropriate wet-dry probability. If the random number is less than or equal to the wet-dry probability, precipitation occurs on that day. Random numbers greater than the wet-dry probability give no precipitation. Since the wet-dry state of the first day is established, the process can be repeated for the next day and so on throughout the simulation period. If wet-dry probabilities are not available, the average monthly number of rainy days may be substituted. The probability of a wet day is calculated directly from the number of wet days: PV = NVD / ND [2.71] where PV is the probability of a wet day, NVD is the number of rainy days, and ND is the number of days, in a month. The probability of a wet day after a dry day can be estimated as a fraction of PV. P(V/D) = ß ?M [2.72] where P(V/D) is the probability of a wet day following a dry day and where /? is a fraction usually in the range of 0.6 to 0.9. The probability of a wet day following a wet day can be calculated directly by using the equation P(V/V) = 1.0 - /? + P(V/D) [2.73] where P(V/V) is the probability of a wet day after a wet day. Vhen /?-4l.O, wet days do not affect probability of rainfall-- p(V/D)=P(V/V)=PV. Conversely, low ß values give strong wet day eifects--/?^0.0, P(V/D)^0., P(V/V)^1.0. Thus, ß controls the interval between rainfall events but has no effect on the number of wet days. For many locations, /?=0.75 gives satisfactory estimates of P(V/D). Although equations 2.72 and 2.73 may give slightly different probabilities than those estimated from rainfall records, they do guarantee correct simulation of the number of rainfall events. 19 Vhen a precipitation event occurs, the amount is generated from a skewed normal daily precipitation distribution: SCF, SCF f (SND. ) (—^) - l)^ - 1-1fin / 6.0 6.0 SCF, RSDVj^ + Rj^ [2.74] where R is the amount of rainfall for day i (mm), SND is the standard normal deviate for day i, SCF is the skew coefficient, RSDV is the standard deviation of daily rainfall (mm), and R is the mean daily rainfall in month k. If the standard deviation and skew coefficient are not available, the model simulates daily rainfall by using a modified exponential distribution. »1= (-In iß) )^\ [2.75] where /¿ is a uniform random number (0.0-1.0) and ( is a parameter usually in the range of 1.0 to 2.0. The larger the ( value, the more extreme the rainfall events. A value oi 1.5 gives satisfactory results at many locations in the united States. The denominator of equation 2.75 assures that the long-term simulated rainfall amount agrees with R. The modified exponential is usually a satisfactory substitute and requires only the monthly mean daily rainfall as input. The amount of daily precipitation is partitioned between rainfall and snowfall according to the average daily air temperature. If the average is O^C or below, the precipitation is snowfall; otherwise, it is rainfall. Air Temperature and Solar Radiation The model developed by Richardson (1981) was selected for use in EPIC because it simulates temperature and radiation, which are mutually correlated with rainfall. The residuals of daily maximum and minimum air temperature and solar radiation are generated from a multivariate normal distribution. The multivariate generation model used implies that the residuals of maximum temperature, minimum temperature, and solar radiation are normally distributed and that the serial correlation of each variable may be described by a first-order linear autoregressive 20 model. Details of the multivariate generation model were described by Richardson (1981). The dependence structure of daily maximum temperature, minimum temperature, and solar radiation was described by Richardson (1982). The temperature model requires monthly means of maximum and minimum temperatures and their standard deviations as inputs. If the standard deviations are not available, the long-term observed extreme monthly minimums and maximums may be substituted. The model estimates standard deviation as 0.25 of the difference between the extreme and the mean for each month. For example, S»™k = »-25 (IE^,k - TiD^.k) P.76] where SDTMX is the standard deviation of the daily maximum temperature, TE is the extreme daily maximum temperature, and T is the average daily maximum temperature for month k. The solar radiation model uses the extreme approach extensively. Thus, only the monthly means of daily solar radiation are required as inputs. The equation for estimating standard deviation is SDRAj^ = 0.25 (RAMXj^ - Rîj^) [2.77] where SDRA is the standard deviation of daily solar radiation (MJ/mi), RAMX is the maximum daily solar radiation at midmonth, and RA is the mean daily solar radiation for month k. Maximum temperature and solar radiation tend to be lower on rainy days. Thus, it is necessary to adjust the mean maximum temperature and solar radiation downward for simulating rainy day conditions. For Tmx this is accomplished by assuming that wet day values are less than dry day values by some fraction of Tmx - TnLmn ' '^V.k = ™«,k - "l (I«,k - ïn,„,k) [2-78] where TV is the daily mean maximum temperature for wet days (oC) in month k, TD is the daily mean maximum temperature for dry days, ftrj, is a scaling factor ranging from 0.0 to 1.0, Tmx is the daily mean maximum temperature, and Tmn is the daily mean minimum temperature. Choosing {lr^=1.0 provides highest deviations on wet days and firp=0.0 ignores the wet day effect. Observed data indicate that firj, usually lies between 0.5 and 1.0. 21 Since equation 2.78 gives lower mean maximum temperature values for wet days, a companion equation is necessary to slightly increase mean maximum temperature for dry days. The development is taken directly from the continuity equation (V,k) ("»k) = (TV,k) (™k) * (™n«,k) ("»»k) P-'»] where ND is the number of days in a month, NVD is the number of wet days, and NDD is the number of dry days. The desired equation is obtained by substituting equation 2.78 into equation 2.79 and solving for TD: NVDk^ . tr . T ^ [2-80] T»mx,k - Tmx,k ^ (^) "T (Tmx,k " T„n,k ) k Use of the continuity equation guarantees that the long-term simulated value for mean maximum temperature agrees with the input value of Tmx- The method of adjusting solar radiation for wet and dry days is similar to that of adjusting maximum temperature. The radiation on wet days is a fraction of the dry day radiation: RAV^ = n^ RADj^ [2.81] where RAV is the daily mean solar radiation on wet days (MJ/m2), iljj is a scaling factor ranging from 0.0 to 1.0, and RAD is the daily mean solar radiation on dry days. An il^ value of 0.5 gives satisfactory results for many locations. The dry day equation is developed by replacing temperature with radiation in equation 2.79 and substituting equation 2.81 for RAV. Then, iU, ) (»D, ) ^^^^^ ^ fij^ (NVDj^ ) + NDDj^ where RA is the daily mean solar radiation for month k (MJ/m2). 22 Vind The wind simulation model was developed by Richardson and Vright (1984) for EPIC. The two wind variables considered are average daily velocity and daily direction. Average daily wind velocity is generated from a two-parameter gamma distribution of the dimensionless form U= (^) 1''^ exp[ (^-1) (1- ^) ] [2.83] Vp Vp where Ü is a dimensionless variable (0-1) expressing frequency with which wind velocity V (m/s) occurs, Vp is the wind velocity at the peak frequency, and TJ is the gamma distribution shape parameter. The shape parameter is calculated with the equation J2_ SDV2 [2.84] where V is the annual average wind velocity (m/s) and SDV is the standard deviation of daily wind velocity (m/s). Values for the average annual wind velocity and the standard deviation of hourly wind are provided by the "Climatic Atlas of the united States" (U.S. Department of Commerce 1968). By experimenting with standard deviations of hourly and daily wind, a correction factor of 0.7 was found to be appropraite for converting hourly standard deviations to daily. The base of the dimensionless gamma distribution (maximum V/Vp) can be determined by Newton's classical method of solving nonlinear equations. The objective function is to select the base to minimize the sum of ln(ü) and 11.5. The value of Vp can be determined by differentiating the gamma function expressed in terms of V and setting the result equal to zero. Then, [2.85] where Vk is the mean daily wind velocity for month k. The rejection technique is used to generate a daily value of V/Vp. The daily wind velocity is then computed by using the equation \ - (V^ ^r^ [2.86] p 23 where Vi is the generated velocity for day i, Vpk is the peak velocity for month k, and V/Vp is the value generated by the rejection technique. Vind direction expressed as radians from north in a clockwise direction is generated from an empirical distribution specific for each location. The empirical distribution is simply the cumulative probability distribution of wind direction. The "Climatic Atlas of the united States" gives monthly percentages of wind from each of 16 directions. Thus, to estimate wind direction for any day, the model draws a uniformly distributed random number and locates its position on the appropriate monthly cumulative probability distribution. Relative Hutidity . 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. The assumed relation between relative humidity on wet and dry days is RHVj^ = MDj^ + flj (1.0 - RHDj^) [2.87] where RHV is the daily mean relative humidity on wet days for month k, RED is the daily mean relative humidity on dry days, and ilg is a scaling factor ranging from 0.0 to 1.0. An llg value of 0.9 seems appropriate for many locations. Using the continuity equation as described in the temperature and radiation sections produces the equation RHk - ng(íí^) RHD, = ÍÍ5_ [2.88] 'k 1.0 - u^f^) ^ ND The where RH is the long-term average relative humidity for month k The appropriate value (RHV or RHD) is used as the peak of a triangular distribution to generate daily relative humidity upper limit of the triangular distribution is set with the equation RHÜ. = RHP. + (1.0 - RHP.) exp(RHP. - 1.0) [2.89] 24 Erosion where RHÜ is the largest relative humidity value that can be generated on day i and RHP is the peak of the triangular distribution (RHV or RED). The lower limit is set with the equation RHL. = RHP. [1.0 - exp(-RHP.) ] [2.90] where RHL is the lowest relative humidity value that can be generated on day i. To assure that the simulated long-term value for mean relative humidity agrees with input RH, the generated value is adjusted by using the equation * RHP. RHG. = RHG. (::^) [2.91] where RHG* is the generated relative humidity on day i adjusted to the mean of the triangle, RHG is. the relative humidity generated from the triangle, and RH is the mean of the triangle. Vater Rainfall/Runoff The EPIC component for water-induced erosion simulates erosion caused by rainfall and runoff and by irrigation (sprinkler and furrow). To simulate rainfall/runoff erosion, EPIC contains three equations--the ÜSLE (Vischmeier and Smith 1978), the MÜSLE iVilliams 1975), 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 depends strictly upon rainfall as an indicator of erosive energy. The MUSLE uses only runoff variables to simulate erosion and sediment yield. Runoff variables increased the prediction accuracy, eliminated the need for a delivery ratio (used in the USLE to estimate sediment yield), and enables the equation to give single storm estimates of sediment yields. 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 uses an equation of the form Y = ;r (K) (CE) (PE) (LS) (ROKF) [2.92] I = El for USLE X = 11.8 (Q* • qp)^-^^ for MUSLE ;t = 0.646 El -f 0.45 (q • q*J^"^^ for Onstad-Foster 25 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, El is the rainfall energy factor, Q* is the runoff volume (m^), q^ is the peak runoff rate (m3/s), Q is the runoff volume (ram), and q*ç is the peak runoff rate (mra/h). The PE value is determined initially by considering the conservation practices to be applied. The value of LS is calculated with the equation (Vischmeier and Smith 1978) LS = (^—) ^ (65.41 S^ + 4.56 S + .065) [2.93] 22.1 where S is the land surface slope (m/m), X is the slope length (m), and ^ is a parameter dependent upon slope. The value of ( varies with slope and is estimated with the equation ^ = 0.3 S / [S + exp(-1.47 - 61.09 S) ] + 0.2 [2.94] The crop management factor is evaluated for all days when runoff occurs by using the equation CE = exp[ (In 0.8 - In CE^^^^.) exp(-1.15 CV) + In CE^^^^^.] [2.95] where CEmn,j 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 0.2 + 0.3 exp(-0.0256 SAN (1 - SIL / 100)1 (1.0 "^^^^ -) (1.0 ( SIL -) 0.3 CLA + SIL 0.25 C C + exp(3.72 - 2.95 C) 0.7 SNl SNl + exp(-5.51 + 22.9 SNl) ■) [2.96] where SAN, SIL, CLA, and C are the sand, silt, clay, and organic carbon contents of the soil (7.) and SN1=1-SAN/100. Equation 2.96 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 26 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>707i). The runoff model supplies estimates of Q and q^. To estimate the daily rainfall energy in the absence of time-distributed rainfall, it is assumed that the rainfall rate is exponentially distributed: ^t " ^p ^^P("t / ^) [2-97] where r is the rainfall rate at time t (mm/h), rp is the peak rainfall rate (mm/h), and k is the decay constant (h). Equation 2.97 contains no assumption about the sequence of rainfall rates (time distribution). The ÜSLE energy equation in metric units is RE = AR (12.1 + 8.9 log —) [2.98] At where RE is the rainfall energy for water erosion equations and AR is a rainfall amount (mm) during a time interval At (h). The energy equation can be expressed analytically as RE = 12.1 /^^ r dt + 8.9 ij" r log r dt [2.99] Substituting equation 2.97 into equation 2.99 and integrating give the equation for estimating daily rainfall energy: RE = R [12.1 + 8.9 (log r - 0.434) ] [2.100] where R is the daily rainfall amount (mm). The rainfall energy factor. El, is obtained by multiplying equation 2.100 by the maximum 0.5-h rainfall intensity (r g) and converting to the proper units: El = R [12.1 + 8.9 (log rp - 0.434) ] (r g) / 1000 [2.101] To compute values for rp, equation 2.97 is integrated to give R = (rp) (/c) [2.102] and R^ = R [1 - exp(-t / K) ] [2.103] 27 The value of R g can be estimated by using a g, as mentioned in the Hydrology section of this chapter: R , = Û , R [2.104] To determine the value of rp, equations 2.104 and 2.102 are substituted into equation 2.103 to give r = -2 R ln(l - o_g) [2-105] Since rainfall rates vary seasonally, « 5 is evaluated for each month by using Veather Service information (U.S. Department of Commerce 1979). The frequency with which the maximum 0.5-h rainfall amount occurs is estimated by using the Hazen plotting position equation (Hazen 1930) p __ \_ [2.106] IT where F is the frequency with which the largest of a total of r 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 Û K, it is necessary to estimate the mean value of R g. The value of R K can be computed easily 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 ^.l,\ .5F,k [2.107] InF, where R . v is the mean maximum 0.5-h rainfall amount, Rgp^j^ is the maximum 0.5-h rainfall amount for frequency F, and subscript k refers to the month. The mean flg is computed with the equation a = ^-^^^ [2.108] 28 where R is the mean amount of rainfall for each event (average monthly rainfall/average number of days of rainfall) and subscript k refers to the month. Daily values of a K are generated from a two-parameter gamma distribution. The base of the gamma distribution is established by examining upper and lower limits of a K- The lower limit determined by a uniform rainfall rate gives a ^ equal to 0.5/24 or 0.0208. The upper limit of Û K 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=l.). The upper limit of a K can be estimated by substituting a high value for rp (250 mm/h is generally near the upper limit of rainfall intensity) into equation 2.103. a 5^ = 1 - exp(-125 / R) [2.109] where a ^ is the upper limit oí a ^. The peak of the a.5 gamma distribution can be computed by using equation 2.85 written in the form ^5,k (^" ^) [2.110] ^5P,k = where 0 ^p 1 is the a K value at the peak of the gamma distribution and 1/ is the gamma distribution shape parameter (a value of 10 is generally satisfactory), and k is the month. The coarse fragment factor is estimated with the equation (Simanton et al. 1984) ROKF = exp(-.03 ROK) [2.111] where ROK is the percent of coarse fragments in the surface soil layer Irrigation Erosion caused by applying irrigation water in furrows is estimated with MUSLE (Villiams 1975): Y = 11.8 (q* • q )^-^^ (K) (CE) (PE) (LS) [2.112] where CE, the crop management factor, has a constant value of 0.5. The volume of runoff is estimated as the product of the irrigation volume applied and the irrigation runoff ratio. 29 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. If irrigation water is applied to land without furrows, the peak runoff rate is assumed to be 0.00189 m3/s per meter of field width. Vind The Manhattan, KS, equation for wind-induced erosion (Woodruff and Siddoway 1965) was modified by Cole et al. (1982) for use in the EPIC model. The original wind erosion equation is of the form VE = f (I, VC, VK, VL, VE) [2.113] where VE is the soil loss from wind erosion (t/ha), I is the soil erodibility index (t/ha), VC is the climatic factor, VK is the soil ridge roughness factor, VL is the field length (mm) along the prevailing wind direction, and VE is the quantity of vegetative cover expressed as small grain equivalent (kg/ha). Equation 2.113 was developed for predicting average annual wind erosion. Its main modification allows EPIC to predict daily values for VE. Two of the variables, I and VC, remain constant for each day of a year. The soil erodibility index is calculated at the start of each year by using a soil textural triangle. Annual I evaluations are necessary to reflect changes in soil surface texture caused by tillage and erosion. The climatic factor, VC, is estimated only once (at the start of a simulation) as described by Lyles (1983). Lyles' method, based on the Thronthwaite precipitation-evaporation index (Thornthwaite 1931) is expressed in the equations VC = 386 V [2.114] Í12 S 10(R - E). i=l and R / 25.4 10/9 10(R - E) = 115 ( ) 1.8 T + 22 [2.115] R > 12.7 mm, T > - I.70C 30 where V is the average annual windspeed (m/s) and R, E, and T are the average values of precipitation (mm), evaporation (mm), and temperature (oC) for month i. The other variables in equation 2.113 are subject to daily variation. The ridge roughness is a function of a row height and row interval 4000 HR^ KR - [2.116] IR where KR is the ridge roughness (mm), HR is the ridge height (m), and IR is the ridge interval (m). The ridge roughness factor is a function of ridge roughness as expressed by the equations VK = 1. , KR < 2.27 ;2.ii7; VK = 1.125 - 0.153 In (KR) , 2.27 < KR < 89. 2.118; VK = 0.336 exp(0.00324 KR) , KR > 89. ;2.ii9; Field length along the prevailing wind direction is calculated by considering the field dimensions and orientation and the wind direction: VL = i^^Um ^2.120] FL I cos(| + G - (j)) I + FV I sin(| + fl - ())) | where FL is the field length (m), FV is the field width (m), 0 is the wind direction clockwise irom north in radians, and (p is the clockwise angle between field length and north in radians. The vegetative cover equivalent factor is simulated daily as a function of the amounts of standing live biomass, standing dead residue, and flat crop residue. VE = 0.2533 (g^ Bj^g + g2 SR + gg FR)^-363 |-2.i2i] 31 where gi, g2, and gs are crop specific coefficients, Bj^g is the above ground biomass of a growing crop (t/ha,) SR is the standing residue from the previous crop (t/ha), and FR is the flat residue (t/ha). Thus, all variables in equation 2.113 can be evaluated. To estimate the soil loss from wind erosion, however, requires a special combination of the factors as follows: E2 = (VK) (I) [2.122] E3 = (VK) (I) (VC) [2.123] VLQ = 1.56 X 10^ (E2)"^-2^ exp(-0.00156 E2) [2.124] VF = E2 [1. - 0.1218 (VL/VLQ)"^-^^2^ exp(-3.33 VL/VL^) [2.125] E4 = (VFÖ-3484 , E3Ö-3484 . E20-3484)2.87 ^2.126] E5 = i>^ E4 ^^2 [2.127] VE = (E5) (DE) / (AE) [2.128] where field lengths greater than VLo (m) do not reduce the erosion estimate, VF is the field length factor, i^i and ^^2 are parameters, DE is the daily wind energy (kVh/m2), and AE is the average annual wind energy (kVh/m2). The parameters rj/i and ^2 are the functions of the vegetative cover factor described by the equations ^2 = 1.+8.93x10"^ VE+8.51x10'^ VE^-1.59x10"^^ VE^ [2.129] i>^ = exp(-7.59x10"^ VE-4.74x10'^ VE^+2.95x10'^^ VE^) [2.130] Daily wind energy is estimated with the equation DE = 193 exp[1.103 C ' ^^ ) ] [2.131] V + 1 where V is the daily average wind velocity (m/s). Average annual wind energy is estimated by integrating the monthly gamma distributions of wind velocity 32 12 AE = 30.4 S k=l (/^ (DE)j^ (;r)k dV ] /^u ;tk dV [2.132] where Vu is the upper limit of wind speed, Vi is the lower limit of erosive wind speed, and x is the frequency of occurrence of wind speed V. Nutrients Nitrogen Nitrate loss in Surface Runoff The amount of NO3-N in runoff is estimated by considering the top soil layer (10-mm thickness) only. The total amount of water leaving the layer is the sum of runoff, lateral subsurface flow, and percolation: qT = Q + 0^ + QR^ [2.133] where QT is the total water lost from the first layer (mm). The amount of NO3-N lost with QT is VN03 = (QT) (cjjgg) [2.134] where VN03 is the amount of NO3-N lost from the first layer and CjTQ« the concentration of NO3-N in the first layer. At the end of the day, the amount of NO3-N left in the layer is VN03 = VN03Q - (QT) (Cj^^g) [2.135] where VN03o and VN03 are the weights of NO3-N contained in the layer at the beginning and ending of the day. The NO3-N concentration can be estimated by dividing the weight of NO3-N by the water storage volume: ^'N03 " ^N03 " ^N03 (^ ^p ^ [2.136] where C.\QO is the concentration of NO3-N at the end of a day, PO is the soil porosity, and VP is the wilting point water content (mm) for soil layer ¿. Equation 2.136 is a finite difference approximation for the exponential equation ^'N03 " ^N03 ^^P(pQ _ yp ) [2.137] 33 Thus, VN03 can be computed for any QT value by integrating equation 2.137: VN03 = VN03 [1 - exp(-^-9T ^ j ^2.138] PO^ - VP^ The average concentration of QT for the day is ^N03 M03 ^2.139] qT Amounts of NO3-N contained in runoff, lateral flow, and percolation are estimated as the products of the volume of water and the concentration from equation 2.139. NO3-N Leaching Leaching and lateral subsurface flow in lower layers are treated by the same approach used for the upper layer except that surface runoff is not considered. NO3-N Transport by Soil Vater Evaporation Vhen water is evaporated from the soil, NO3-N is moved upward into the top soil layer by mass flow. The equation for estimating this NO3-N transport is EN03 = S SEV¡ (CjjQg)^ [2.140] where EN03 is the amount of NO3-N (kg/ha) moved from lower layers to the top layer by soil water evaporation Eg (mm), subscript ¿ refers to soil layers, and M is the number of layers contributing to soil water evaporation (maximum depth is 0.2 m). Organic N Transport by Sediment A loading function developed by McElroy et al. (1976) and modified by Williams and Hann (1978) for application to individual runoff events is used to estimate organic N loss. The loading function is YON = 0.001 (Y) (CQJJ) (ER) [2.141] where YON is the organic N runoff loss (kg/ha), Y is the sediment yield (t/ha), Cr.^ is the concentration of organic N in the top soil layer (g/t), and ER is the enrichment ratio. The enrichment ratio is the concentration of organic N in the sediment divided 34 by that in the soil. Enrichment ratios are logarithmically related to sediment concentration as described by Menzel (1980). An individual event enrichment-sediment concentration relationship was developed for EPIC considering upper and lower bounds. The upper bound of enrichment ratio is the inverse of the sediment delivery ratio. Exceeding the inverse of the delivery ratio implies that more organic N leaves the watershed than is dislodged from the soil. The delivery ratio is esti- mated for each runoff event by using the equation q„ 0,56 DR = (-fi-) [2.142] where DR is the sediment delivery ratio (sediment yield divided by gross sheet erosion), q*p is the peak runoff rate (mm/h), and rep is the peak rainfall excess rate (mm/h). Equation 2.142 is based on sediment yield estimated by using MÜSLE (Villiams 1975). The rainfall excess rate cannot be evaluated directly because tne hydrology model predicts only the total daily runoff volume. An estimate of the rate can be obtained, however, using the equation ^ep = ^p - f C2.143] where rp is the peak rainfall rate (mm/h) and f is the average infiltration rate (mm/h). The average infiltration rate can be computed from the equation f = ^-^ [2.144] DUR where DUR is the rainfall duration (h). The rainfall duration can be estimated by solving equations 2.102 and 2.103 for t when Rt/R=0.99. Thus, DUR = 4^605_L [2.145] r P The lower limit of enrichment ratio is 1.0--sediment particle size distribution is the same as that of the soil. Thus, 1< ER< 1/DR. The logarithmic equation for estimating enrichment ratio is x^ ER = x^ Cg "^ [2.146] 35 where Cs is the sediment concentration (g/m^) and xi and X2 are parameters set by the upper and lower limits. For the enrichment ratio to approach 1.0, the sediment concentration must be extremely high. Conversely, for the enrichment ratio to approach 1/DR, the sediment concentration must be low. The simultaneous solution of equation 2.146 at the boundaries assuming that sediment concentrations range from 500 to 250,000 g/m^ gives X. = -log(^ / 2.699 [2.147] ^ DR Xw = ^ (0.25)''2 [2.148] Dentrification As one of the microbial processes, denitrification is a function of temperature and water content. The equation used to estimate the denitrification rate is DN¿ = VN03^ [l - exp[(-1.4) (TFj^^ (C^ ]j, SVF > 0.9 [2.149] DN = 0. , SVF < 0.9 DN is the denitrification rate in layer I (kg/ha/d), TFn is the nutrient cycling temperature factor, C is the organic carbon content (7.), and SVF is the soil water factor. The temperature factor is expressed by the equation ^ , T^ > 0. [2.150] ^\t T^ + exp(9.93 - 0.312 T^) TFjj^ = 0., T^ < 0- where T is soil temperature (oC) and subscript I refers to the layers. The soil water factor considers total soil water in the equation SV. SVF, = —^ [2.151] ^ PO^ where SV is the soil water content in layer I and PO is the soil porosity (mm). 35 lineralization The N mineralization model is a modification of the PAPRAN mineralization model (Seligman and van Keulen 1981). The model considers two sources of mineralization: fresh organic N pool, associated with crop residue and microbial biomass, and the stable organic N pool, associated with the soil humus. Mineralization from the fresh organic N pool is estimated with the equation RMN^ = (DCR^) (FON^) [2.152] where RMN is the N mineralization rate (kg/ha/d) for fresh organic N in layer ¿, DCR is the decay rate constant for the fresh organic N, and FON is the amount of fresh organic N present (kg/ha). The decay rate constant is a function of C:N ratio, C:P ratio, composition of crop residue, temperature, and soil water: DCR. = (CNP.) (RC) (—^) • TF^.)^-^ [2.153] where CNP is a C:N and C:P ratio factor, RC is a residue composition factor, and FC is the soil water content (mm) at field capacity. The value of CNP is calculated with the equation (exp[-0.693 (CNR - 25) / 25] CNP, . min exp[-0.693 (CPR^ - 200) / 200] ^2.154] ^ ll.O where CNR is the C:N ratio and CPR is the C:P ratio in layer ¿. The C:N and C:P ratios of crop residue are computed for each soil layer with the equations 0.58 FR. CNR, = ^ [2.155] ^ FON^ + VN03^ 0.58 FR, CPR, = ^ [2.156] FOP^ + AP^ where FOP is the amount of fresh organic P in layer ¿ (kg/ha) and AP is the amount of labile P (kg/haj. The value of RC is determined by the stage of residue decomposition. The first 20X of the residue is decomposed by using RC=0.8 (a rate appropriate for carbohydrate-like material). Between 20 and 907i, RC=0.05 (cellulose-like material). The final 107. of the residue is decomposed at a rate appropriate for lignin (RC=0.0095). 37 Organic N associated with humus is divided into two pools--active and stable--by using the equation ON^^ = (RTN^) (ON^) [2.157] where ONa is the active or readily mineralizable pool (kg/ha), RTN is the active pool fraction, ON is the total organic N (kg/ha), and the subscript I is the soil layer number. The active pool fraction in the plow layer depends on the number of years the soil has been cultivated and is estimated with the equation RTN^ = 0.4 exp(-0.0277 YC) + 0.1 [2.158] where YC is period of cultivation before the simulation starts (yr). The concepts expressed in equation 2.158 are based on work of Hobbs and Thompson (1971). Below the plow layer the active pool fraction is set to 407i of the plow layer value, based on work of Cassman and Munns (1980). Organic N flux between the active and stable pools is governed by the equilibrium equation RON, = BKN ON . (^—) - ONÍ = ^^^^ r^i ^—f - "^¿j [2.159] where RON is the flow rate (kg/ha/d) between the active and stable organic N pools, BKN is the rate constant («1 X 10"5) (d"0, ON is the stable organic N pool, and subscript P is thes soil layer number. The daily flow of humus related organic N (RON) is added to the stable pool and subtracted from the active pool. Only the active pool of organic N is subjected to mineralization. The humus mineralization equation is HMN^ = (CMN) (ON^p (SVF^ • TFj^^)^-^ (BD^)^ / (BDP^)^ [2.160] where HMN is the mineralization rate (kg/ha/d) for the active organic N pool in layer ¿, CMN is the humus rate constant «0.0003) (d-i), BD is the settled bulk density of the soil t/m3), ana BDP is the current bulk density as affected by tillage (t/m3). To maintain the N balance at the end of a day, the humus mineralization is subtracted from the active organic N pool; the residue mineralization is subtracted from the FON pool; 207. of RMN is added to the active ON pool; and 807. of RMN is added to VN03 pool. 38 Immobilization Like the mineralization model, the immobilization model is a modification of the PAPRAN model. Immobilization is a very important process in EPIC because it determines the residue decomposition rate. Of course, residue decomposition has an important effect on erosion. The daily amount of immobilization is computed by subtracting the amount of N contained in the crop residue from the amount assimilated by the microorganisms: VIM^ = (DCR^) (FR^) (0.016 - c^p^) [2.161] where VIM is the N immobilization rate in layer Í (kg/ha/d); 0.016 is the result of assuming that C=0.4 FR, that C:N of the microbial biomass and their labile products = 10, and that 0.4 of C in the residue is assimilated; and Cjrpp is the N concentration in the crop residue (g/g). Immobilization may be limited by N or P availability. If the amount of N available is less than the amount of immobilization predicted from equation 2.161, the decay rate constant is adjusted with the relationship 0.95 VN03. DCR'. = [2.162] FR^ (0.016 - Cjjpj) where DCR' allows 957. use of the available NO3-N in soil layer Í. A similar adjustment is made if P is limiting. The crop residue is reduced by using the equation FR^ = FR^^ - (DCR'^) (FR^^) [2.163] where FRQ and FR are the amounts of residue in soil layer i at the start and end of a day (kg/ha). Finally, the immobilized N is added to the FON pool and subtracted from the VN03 pool. Rainfall To estimate the N contribution from rainfall, EPIC uses an average rainfall N concentration for a location for all storms. The amount of N in rainfall is estimated as the product of rainfall amount and concentration. Phosphorus Soluble P loss in Surface Runoff The EPIC approach is based on the concept of partitioning pesticides into the solution and sediment phases as described by Leonard and Vauchope (Knisel 1980). Because P is mostly associated with the sediment phase, the soluble P runoff equation can be expressed in the simple form 39 YSP = 0.01 (Cj^p^) (Q) / k^ [2.164] where YSP is the soluble F (kg/ha) lost in runoff volume Q (mm), Crp/, is the concentration oi AP in soil layer Í (g/t), and kd is the F concentration of the sediment divided by that of the water (m3/t). The value of kd used in EFIC is 175. P Transport by Sediment Sediment transport of F is simulated with a loading function as described in organic N transport. The F loading function is YF = 0.001 (Y) (c ) (ER) [2.165] where YF is the sediment phase F lost in runoff (kg/ha) and Cp is the concentration of F in the top soil layer (g/t). lineralization The F mineralization model developed by Jones et al. (1984) is similar in structure to the N mineralization model. Mineralization from the fresh organic F pool is estimated for each soil layer with the equation RMF^ = (DCR^) (FOF^) [2.166] where RMF is the mineralization rate of fresh organic F in layer I (kg/ha/d) and FOF is the fresh organic F in crop residue (kg/ha). Mineralization of organic F associated with humus is estimated for each soil layer by using the equation „P ("«») i^\l) (^h) (SVf£ ♦ ^hl)"-" (BPP' p^^^^j where HMF is the humus F mineralization rate (kg/ha/d) and OF is the organic F content of soil layer i (kg/ha). The rate of ONa conversion to ON is used in equation 2.167 to calculate the active portion of the OF pool. This eliminates the need for maintaining two OF pools corresponding to the ONa and ONg pools. To maintain the F balance at the end of a day, humus mineralization is subtracted from the organic F pool; residue mineralization is subtracted from the FOF pool; 207i of RMF is added to the OF pool; and 807. of RMF is added to labile F pool. 40 Immobilization The P immobilization model, also developed by Jones et al. (1984), is similar in structure to the N immobilization model. The daily amount of immobilization is computed by subtracting the amount of P contained in the crop residue from the amount assimilated by the microorganisms: VIP^ = (DCR'^) (FR^) (0.16 LFj^ - Cp^^) [2.168] where VIP is the P immobilization rate in layer t (kg/ha/d); 0.16 is the result of assuming that C=0.4 FR and that 0.4 of the C in FR is assimilated by soil microorganisms; LFj is the labile P immobilization factor allowing the P:C ratio of soil microorganisms to range from 0.01 to 0.02 as a function of labile P concentration, and Cppi. is the P concentration in the crop residue. The labile P immobilization factor is computed with the equations LF U 0.01 + 0.001 c LP¿ \^l ^ 10 LFj^ =0.02 ^LP¿ > 10 [2.169] [2.170] The immobilized P is added to the FOP pool and subtracted from the labile P pool. lineral F Cycling The mineral P model was developed by Jones et al. (1984). Mineral P is transferred among three pools: labile, active mineral, and stable mineral. Fertilizer P is labile (available for plant use) at application but may be quickly transferred to the active mineral pool. Flow between the labile and active mineral pools is governed by the equilibrium equation PSP MPR^ = 0.1 SVF^ exp(0.115 T^ - 2.88) [AP^ - MP^^ ( ^)| [2.171] where MPR is the mineral P flow rate for layer i (kg/ha/d), T is the soil temperature (^C), MPa is the amount in the active mineral P pool fkg/ha), and PSP is the P sorption coefficient defined as the traction of fertilizer P remaining in the labile pool after the initial rapid phase of P sorption is complete. The daily amount of P computed with equation 2.171 flows to the active mineral P pool and is, therefore, added to that pool and subtracted from the labile pool. Obviously, the flow reverses when labile P is less than MP^ PSP^ / (1-PSP^). The P sorption 41 coefficient is a function of chemical and physical soil properties as described by the following equations (Jones et al. 1984). In calcareous soils PSP^ = 0.58 - 0.0061 CAC^ [2.172] In noncalcareous, slightly weathered soils PSP^ = 0.02 + 0.0104 AP^ [2.173] In noncalcareous, moderately weathered soils PSP^ = 0.0054 BSA^ + 0.116 PH^ - 0.73 [2.174] In noncalcareous, highly weathered soils PSP^ = 0.46 - 0.0916 In CLA^ [2.175] where PSP is the P sorption coefficient for soil layer ¿, CAC is the CaCOa concentration (g/t^, and BSA is the base saturation by the ammonium acetate (NH4OACJ method (X). PSP is constrained within the limits 0.05<PSP<0.75. At equilibrium the stable P pool is assumed to be four times as large as the active mineral P pool. Flow between the P pools is governed by the equation ASPR^ = i^¿ (4 MP^^ - MPg^) [2.176] where ASPR is the flow rate between the active and stable mineral P pools (kg/ha/d) for soil layer ¿, w is the flow coefficient 'd-i), and IPs is the amount of stable mineral P (kg/ha). The ally amount of P computed with equation 2.176 flows into the stable pool and is subtracted from the active pool. Obviously, the flow reverses when MP^>4MP^. The flow coefficient, á;, is a function of PSP as expressed by the equations (Jones et al. 1984) ï tí^ = exp(-1.77 PSP^ - 7.05) [2.177] for noncalcareous soils, and uj¿ = 0.0076 [2.178] for calcareous soils. 42 Soil Temperature Daily average soil temperature at the center of each soil layer is simulated for use in nutrient cycling and hydrology. The basic soil temperature equation is T^^. = LAG(T^^._^) + (1.0 - LAG) (Z^ (T ^ TG.) + TG.) [2.179] where T is the soil temperature at the center of layer I on day i (oC), LAG is a coefficient ranging from 0.0 to 1.0 that allows proper weighting of yesterday's temperature, T is the long-term average annual air temperature at the site, TG is the soil surface temperature, and FZ is a depth factor. Thus, given yesterday's temperature, equation 2.179 estimates today's temperature as a function of soil surface temperature, depth, and a lag coefficient. It is assumed that the temperature remains almost constant at some depth called damping depth and is approximately T. The depth weighting factor governs temperature changes between the soil surface and the damping depth according to the equation FZ. = [2.180] ^ ZD + exp(-0.867 - 2.08 ZD) where Z^ + Z^ ^ ZD = -^ ^^^^ [2.181] 2.0 DD where Z is soil depth from the surface (m) and DD is the damping depth (m). Obviously, equations 2.180 and 2.181 make near surface temperatures a strong function of TG. As depth increases, T has more influence until finally at the damping depth, the temperature is within 57i of T. The damping depth is a function of soil bulk density and water content as expressed in the equation DP = 1.0 + 2^^_M [2.182] § = — [2.183] BD + exp(6.53 - 5.63 BD) SV (0.356 - 0.144 BD) Zjj 2 DD = DP exp iln(^) ir^-^ \ [2.184] ^ DP 1 + § ^ 43 where DP is the maximum damping depth for the soil (m), BD is the soil bulk density (t/m^), Zj, is the distance from the bottom of the lowest soil layer to the surface, and § is a scaling parameter. To complete the solution of equation 2.179, the soil surface temperature must be estimated. The first step is to estimate the bare soil surface temperature (TGB). Of course, TGB is usually closely related to the air temperature. Other important factors that also influence TGB are precipitation and previous soil temperature. Vhen precipitation occurs, the soil surface temperature usually decreases. Thus, the appropriate air temperature for estimating TGB is near the daily minimum. TGBV. = T . + n (T . - T .) [2.185]1 mn,i s ^ mx,i mn,i^ ^ -^ where TGBV is the bare soil surface temperature on wet day, i, and ils is a scaling factor to adjust for wet days. The value of fig ranges from 0.0 to 1.0, but more realistic results can be obtained using fis«0.1. The companion equation for dry days is derived from the continuity equation (equation 2.79). T . + T NVD. ^^^ ^^^^ _ ( íi^ TT + ft ÍT - T )^ o Nn "^^ji s "^ji mn,i^-J TGBD. "»k NVD, 1.0 «»k [2.186] where TGBD is the bare soil surface temperature on dry day, i, NVD is the number of wet days, and ND is the number of days in month k. To estimate the lag in the system caused by heat stored in the soil, a 5-day moving average is applied to TGB. * 4 TGB. = S TGB. [2.187] ^ N=0 '"^ where TGB* is the final estimate of bare soil surface temperature (oC) and TGB is either TGBV or TGBD obtained from equations 2.185 and 2.186. 44 If the soil surface is not bare, the surface temperature can be affected considerably by the amount of cover (crop residue or snow). This effect can be simulated by lagging the predicted bare surface temperature according to the equation TG. = (bcv) (TGB*._^) + (1 - bcv) (TGB*.) [2.188] where TG is the final estimate of soil surface temperature (oC) and bcv is a lagging factor for simulating residue and snow cover effects on surface temperature. The value of bcv is 0 for bare soil and approaches 1.0 as cover increases, as expressed in the equation CVr CV + exp(7.563 bcv = max I CV + exp(7.563 - 1.297 X 10-4 CV) [2.189] SNO SNO + exp(2.303 - 0.2197 SNü) where CV is the sum of above ground biomass and crop residue (t/ha) and SNO is the water content of the snow cover (mm). Crop Growth Model A single model is used in EPIC for simulating all the crops considered (Corn, Grain sorghum, Wheat, Barley, Oats, Sunflower, Soybean, Altalfa, 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. EPIC is capable of simulating growth for both annual and perennial crops. 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 they may become dormant after frost. They start growing when the average daily air temperature exceeds their base temperature. Phenological development of the crop is based on daily heat unit accumulation. It is computed by using the equation T + T HUj^ = íjnxik mA) _ T^ . HUj^ > 0 [2.190] 45 where Hü, T^x, and Tnm are the values of heat units, maximum temperature, and minimum temperature (»C) on day k, and Tb is the crop-specific base temperature (oC) (no growth occurs at or below Tb) of crop j. A heat unit index (HÜI) ranging from 0 at planting to 1 at physiological maturity is computed as follows: I i ^ E k=l HÜ, HÜI. = ^ PHÜ. [2.191] where HÜI is the heat unit index for day i and PHÜ is the potential heat units required for the maturation of crop j. The value of PHÜ may be inputted or calculated by the model from normal planting at harvest dates. Date of harvest, leaf area growth and senescence, optimum plant nutrient concentrations, and partition of dry matter among roots, shoots, and economic yield are affected by HÜI. Potential Growth Interception of solar radiation is estimated with a Beer's law equation (Monsi and Saeki 1953) PAR. =0.5 (RA). [1. - exp(-0.65 LAI) ]. [2.192] where PAR is photosynthetic active radiation (MJ/m^), RA is solar radiation (MJ/m^), LAI is the leaf area index, and subscript i is the day of the year. Using Monteith's approach (Monteith 1977), potential increase in biomass for a day can be estimated with the equation AB .= 0.001 (BE). (PAR). (1 + AHRLT.)^ [2.193] where ABp is the daily potential increase in biomass (t/ha), BE is the crop parameter for converting energy to biomass ikg/MJ), HRLT is the day length (h), and AHRLT is the change in day length (h/d). The day length tunction of equation 2.193 increases potential growth during the spring and decreases it in the fall (Baker et al. 1980). Day length is a function of the time of year and latitude as expressed in the equation HRLT. = 7.64 cos~^ -tan(^LAT) tan(SD). [2.194] ^ ^ 365 ^^ 45 where LAT is the latitude of the watershed (degrees) and SD, the sun's declination angle, is defined by the equation 2T SD. = 0.4102 sin[ ^^ (i ^ 365 80.25) ] [2.195] In most crops, leaf area index (LAI) is initially zero or very small. It increases exponentially during early vegetative growth, when the rates of leaf primordia development, leaf tip appearance, and blade expansion are linear functions of heat unit accumulation (Tollenaar et al. 1979; Vatts 1972). In vegetative crops such as sugarcane and some forages, LAI reaches a plateau, at which time the rates of senescence and growth of leaf area are approximately equal. In many crops, LAI decreases after reaching a maximum and approaches zero at physiological maturity. In addition, leaf expansion, final LAI, and leaf duration are reduced by stresses (Acevedo et al. 1971; Eik and Hanway 1965). LAI is simulated as a function of heat units, crop stress, and crop development stages. From emergence to the start of leaf decline, LAI is estimated with the equations LAI. = LAI.^ + ALAI [2.196] ALAI = (AHUF) (LAI^) 1. - exp[5.0(LAI._^ - Ul^) ] y REG. . [2.197] where LAI is leaf area index, HÜF is the heat unit factor, and REG is the value of the minimum crop stress factor. Subscript mx is the maximum value and A is the daily change. The heat unit factor is computed by using the equation HÜF. = HÜI. HUI . + exp[ah.^^ - (ah. ^2) (^'^^i) ^ [2.198] where ahj ,1 and ahj ,2 are parameters of crop j, and HÜI is the heat unit index. From the start of leaf decline to the end of the growing season, LAI is estimated with the equation LAI. = LAI^ Í1 - HÜI.1 1 - HÜI, ad. [2.199] 47 where ad is a parameter that governs LAI decline rate for crop j and subscript o is the day of the year when LAI starts declining. Crop height is estimated with the relationship CHT. = HMX. /ÎUFT [2.200] where CHT is the crop height (m) and HMX is the maximum height for crop j. 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 simulates this partitioning by decreasing tne fraction linearly from 0.4 at emergence to 0.2 at maturity. Thus, the potential daily change in root weight is computed with the equation ARVT. = ABp . (0.4 - 0.2 HÜI.) [2.201]JL 1 • X JL where ARVT is the change in root weight (t/ha) on day i. The potential change in root weight through the root zone is simulated as a function of plant water use in each layer of soil with the equation Í ^ J^ 1 ARV.^^ = (ARVT.) Uj [2.202]H E u. f where RV is the root weight in soil layer Í (t/ha), M is the total number of soil layers, and u is the daily water use rate in layer ¿ (mm/d). Rooting depth normally increases rapidly from the seeding depth to a crop-specific maximum. In many crops, the maximum is usually attained well before physiological maturity (Borg and Grimes 1986). Rooting depth is simulated as a function of heat units and potential root zone depth: ARD. =2.5 (RDMX.) (AHUF.) , RD. < RZ. [2.203] where RD is the root depth (m), RDMX is the maximum root depth 'm) for crop j in ideal soil, and RZ is the soil profile depth m). The economic yield of most grain, pulse, and tuber crops is a reproductive organ. Crops have a variety of mechanisms which ensure that their production is neither too great to be supported 43 by the vegetative components nor too small to ensure survival of the species. As a result, harvest index (economic yield/above- ground biomass) is often a relatively stable value across a range of environmental conditions. In EPIC, crop yield is estimated by using the harvest index concept: YLD. = (HIj) (B^Q) [2.204] where YLD is the amount of the crop removed from the field (t/ha), HI is the harvest index, and BJ^Q is the above-ground biomass (t/ha) for crop j. Harvest index increases nonlinearly from 0 at planting to 1.0 at maturity according to the equation HIA. = HI. i ^ S AHUFH, k=l ^ [2.205] where HIA is the harvest index on day i and HUFH is the heat unit factor that affects harvest index. The harvest index heat unit is computed with the equation HÜI. HUFH. = [2.206] ^ HÜI. + exp(6.50 - 10.0 HUI^) The constants in equation 2.206 are set to allow HUFHi to increase from 0.1 at HÜIi=0.5 to 0.92 at HÜIi=0.9. This is consistent with the economic yield development of grain crops, which produce the greatest economic yield in the second half of the growing season. Vater Use The potential water use, Ep, is estimated as a fraction of the potential evaporation by using the leaf-area-index relationship developed by Ritchie (1972). LAI. Epi = \i (^) ' hi <- \i t2.207] where EQ is the potential evaporation and LAI is the leaf area index on day i. 49 The potential water use from the soil surface to any root depth is estimated with the function E Pi ^Pi exp(-A) 1. - exp[-A (^)] RZ [2.208] where Up is the total water use rate (mm/d) to depth Z (m) on day i, RZ is the root zone depth (m), and A is a water use distribution parameter. The amount used in a particular layer can be calculated by taking the difference between Up^i^ values at the layer boundaries: 'Pi '?¿ exp(-A) exp[-A (-^)] RZ 1. exp[-A (^] RZ [2.209] where Upy? is the potential water use rate from layer ¿ (mm/d). Equation 2.209 applies to a soil that provides poor conditions for root development when A is set to a high value like 10. The high A value gives high water use near the surface and very low use in the lower half of the root zone. Since there is no provision for water deficiency compensation in any layer, considerable water stress may be incorrectly indicated if equation 2.209 is used. To overcome this problem, equation 2.209 was modified to allow plants to compensate for water deficiency in a layer by using more water from other layers. Total compensation can be accomplished by taking the difference between Up. at the bottom of a layer and the sum of water use above a layer: ^Pi '?¿ exp(-A) exp[-A (^ ) ] RZ M S uj K=l * [2.210] where Uk is the actual water use rate (mm/d) for all layers above layer ¿. Thus, any deficit can be overcome if a layer that is encountered has adequate water storage. Neither equation 2.209 (no compensation) nor equation 2.210 (total compensation) is satisfactory to simulate a wide range of soil conditions. A combination of the two equations, however, provides a very general water use function: %i '?£ exp(-A) 50 1.- exp[-A (-i)] - (1. - ÜC) RZ 1. - exp[-A (-^)] RZ ^ í-1 - ÜC S u, k=l * [2.211] where ÜC varies over a range (O.-l.) and is the water deficit compensation factor. In soils with a good rooting environment, ÜC=1. gives total compensation. The other extreme, poor conditions, allow no compensation (ÜC=0.). The procedure for estimating ÜC is described in the Growth Constraints section of this chapter. The potential water use in each layer calculated with equation 2.211 is reduced when the soil water storage is less than 257. of plant-available soil water (Jones and Kiniry 1986) by using the equation 'YÍ exp 4. (SV.. - VP.) 5.1 — — - 1. (FC¿ - VPp FC. - VF SV^< ¿ VF í 'f¿ FC, - VF. sv¿ > —^ — + VP^ [2.212] [2.213] where SV is the soil water content in layer ¿ on day i (mm) and FC and VF are the soil water contents at field capacity and wilting point for layer Í. Nutrient Uptake Nitrogen Supply and Demand. Crop use of N is estimated by using a supply and demand approach. The daily crop N demand is the difference between the crop N content and the ideal N content for that day. The demand is estimated with the equation ÜNDi = (CNB)í (B)i - V ÜNj^ [2.214] 51 where UND is the N demand rate of the crop (kg/ha/d), Cj^g is the optimal N concentration of the crop (kg/t), B is the accumulated biomass (t/ha] for day i, and UN is the actual N uptake rate (kg/ha/d). Tne optimal crop N concentration declines with increasing growth stage (Jones 1983a) and is computed as a function of growth stage by using the equation ^NBi " ^^1 "*" ^^2 ^^P("^^3 ^^^i) [2.215] where bni, bn2, and bna are crop parameters expressing N concentration and HÜI (heat unit index) is the fraction of the growing season. Soil supply of N is assumed to be limited by mass flow of NO3-N to the roots ÜN^,i - ^t,i iVN03^1 SV^ [2.216] where UN is the rate of N supplied by the soil (kg/ha/d), VN03 is the amount of NO3-N (kg/ha), SV is the soil water content (mm), u is water use rate (mm/d), and subscript i refers to the soil layers. The total mass flow demand is estimated by summing the layer demands: H UNS, = S UN. . [2.217] where UNS is the N supply rate from soil to plants (kg/ha/d). Since mass flow uptake can produce questionable results when N concentrations are extremely high or low, UN values obtained from equation 2.216 are adjusted: UND. ÜNa. . = UN. . ( -) , ÜNa. . < VN03. • [2.218] Equation 2.218 assures that actual N uptake cannot exceed the plant demand when mass flow estimates are too large. It also provides for increased N supply when mass flow estimates are too low despite the availability of NO3. Fixation. Daily N fixation is estimated as a fraction of daily plant N uptake for legumes: VFX. = FXR. • UN. , VFX < 6.0 [2.219] 52 where VFX is the amount of N fixation (kg/ha) and FXR is the fraction of uptake for day i. The fraction, FXR, is estimated as a function of soil NO3 and soil water contents and plant growth stage: FXR = min(1.0, FXV, FXN) • FXG [2.220] where FXG is the growth stage factor, FXV is the soil water content factor, and FXN is the soil NO3 content factor. The growth stage factor is computed with the equations FXG. = 0.0 , HÜI. < 0.15 , HUI^ > 0.75 [2.221] FXG. = 6.67 HÜI. - 1.0 0.15 < HUI. < 0.3 [2.222] FXG. = 1.0 0.3 < HUI. < 0.55 [2.223] FXG. = 3.75 - 5.0 HUI. 0.55 < HUI. < 0.75 [2.224] where HUI is the heat unit index for day i. The soil water content factor reduces N fixation when the water content in the top 0.3 m is less than 857. of field capacity according to the equation SV3. - VP3 FXV. = i , SV3 < 0.85(FC3 - VP3) + VP3 [2.225] ^ 0.85 (FC3 - VP3) where SV3, VP3, and FC3 are the water contents in the top 0.3 m of soil on day i, at wilting point, and at field capacity. The amount of NO3 in the root zone determines the soil NO3 factor, FXN: FXN = 0. , VN03 > 300. kg/ha/m [2.226] FXN = 1.5 - 0.005 {^^^) , 100 < VN03 < 300. [2.227] RD FXN = 1.0 , VN03 < 100. kg/ha/m [2.228] where VN03 is the weight of N03 in the root zone (kg/ha) and RD is the root depth (m). 53 Phosphorus Crop use of P is estimated by the supply and demand approach described in the N model. The daily plant demand is computed with equation 2.214 written in the form üPDi = (CpB)i(B)i - 3 üPj^ [2-229] where ÜPD is the P demand rate for the plant (kg/ha/d), Cpg is the optimal P concentration for the plant, and UP is the actual P uptake rate (kg/ha/d). The optimal plant P concentration is computed with equation 2.215 written in the form Cpg. = bp^ ^ bp2 exp(-bp3 HÜI.) [2.230] where bpi, bp2 5 and bp3 are crop parameters expressing P concentration. Soil supply of P is estimated by using the equation M RV UPSj = 1.5 UPDj ^S^ (IPJ^ (-p|^) [2.231] where UPS is the rate of P supplied by the soil (kg/ha/d), LFu is the labile P factor for uptake, RV is the root weight in layer Í (t/ha), and RVT is the total root weight on day i (t/ha). The constant 1.5 allows two-thirds of the roots to meet the P demand of the plant if labile P is not limiting. The labile P factor for uptake ranges from 0.1 to 1.0 according to the equation ^^Mi = ^-^ ^ Cj^p^ -f 117. exp(-0.283 Cj^p^) ^2-232] where Crp is the labile P concentration in soil layer Í (g/t). Equation 2.232 allows optimum uptake rates when Cj^p is higher than 20 g/t. This is consistent with critical labile P concentrations for a range of crops and soils (See Chapter 7). Sharpley et al. (1984, 1985) described methods of estimating Cj^p from soil test P and other soil characteristics. Growth Constraints usually the potential crop growth and yield are not achieved because of various constraints imposed by the plant environment. 54 The model estimates the severities of the stresses caused by water, nutrients, temperature, aeration, and radiation. The estimates (stress factors) range from 0.0 (the most severe) to 1.0, and the stresses affect plants in several ways. In EPIC, the stresses are considered in estimating constraints on biomass accumulation, root growth, and yield. The biomass constraint is calculated by using the lowest value from among the stress factors estimated for water, nutrients (N and P), temperature, and aeration. The root growth constraint is the minimum of the soil strength, temperature, and aluminum toxicity stresses. A description of the stress factors involved in determining each constraint follows. Biomass The potential biomass predicted with equation 2.193 is adjusted daily with the following equation if any one of five plant stress factors is less than 1.0: AB = (ABp) (REG) [2.233] where REG is the crop growth regulating factor (the lowest value from among the estimates for the stress factors). Vater Stress. The water stress factor is computed by considering supply and demand in the equation H Su., VSi - -¥^ [2.234] Pi where VS is the water stress factor, u is the water use in layer ^, and Ep is the potential plant water evaporation rate on day i. This is consistent with the concept that drought stress limits biomass production in proportion to transpiration reduction (Hanks 1983). Temperature Stress. The plant temperature stress factor is estimated with the equation TS. ^ sin U (T , . g, ) [2.235] ^ oj bj ^ where TS is the plant temperature stress factor, TG is the soil surface temperature (oC), Tb is the base temperature for crop j, and To is the optimal temperature for crop j. Equation 2.235 produces symmetrical plant growth stress about the optimal temperature and considers average daily soil surface temperature as a stress indicator. 55 Kutrient Stress. The N and P stress factors are based on the ratio of accumulated plant N and P to the optimal values. The stress factors vary nonlinearly from 1.0 at optimal N and P levels to 0. when N or P is half the optimal level (Jones 1983a) For N, the scaling equation is '\i - 2 S UN, k=l * í^NB)í ^»)i [2.236] where SNg is a scaling factor for the N stress factor, UN is the crop N uptake rate on day k (kg/ha/d), Cj^g is the optimal N concentration of the crop on day i, and B is the accumulated biomass (t/ha). The N stress factor is computed with the equation SN. = 1 ^^^ [2.237] ^ SNg . + exp(3.39 - 10.93 SNg •) where SN is the N stress facto