The effect of soil and plant water deficits on plant growth and yield reduction is calculated by the soil water balance. The soil water balance routine of CERES-Wheat includes the soil water quantity resulting from the input of precipitation and irrigation, the outputs of evaporation from plants and the soil, and runoff and drainage. In the model, the soil water balance is distributed in up to ten layers. The user assigns the depth of each layer. The water balance routines described in this chapter were not designed to be accurate for wetlands where a water table is present near the surface nor where artificial drainage occurs. It was primarily developed to evaluate naturally drained soils where shortages in the soil water balance are likely to occur, not where aeration problems may reduce crop growth.
The water content of any particular soil layer can decrease as a result of soil evaporation, root absorption, or water flow to an adjacent layer. The limits to which water can increase or decrease are inputs for each soil layer. These limits are the lower limit of plant water availability (LL(L)), the field drained upper limit (DUL(L)), and the field saturated water content (SAT(L)). When the crop water supply is marginal, the accuracy of these input limits becomes quite important for accurate prediction of crop performance.
Traditional laboratory measurements of wilting point and field capacity water content are frequently inaccurate for establishing field limits of water availability (Ritchie 1981). Field measurements are preferred when possible. A comparison of field-measured versus laboratory-measured limits of soil water content values taken from several sites in the United States revealed that considerable error was possible when laboratory methods were applied to the field limits (Ratliff et al. 1983).
Other soil inputs needed in the model include the soil albedo (SALB), the upper limit of first stage soil evaporation (U), a constant for calculating the drainage rate (SWCON), and a curve number used to calculate run-off (CN2).
If irrigation water is applied to the field being simulated, then the day of the year of irrigation (JDAY(J)) and the amount of irrigation (AIRR(J)) are read for the number of irrigations (NIRR) specified in the data. When the soil water balance routines begin on the date specified by the beginning of weather data, several variables have to be initialized with some reasonable logic. This initialization is done in the subroutine where soil properties are read from a file and initialized (SOILRI).
INITIATION OF PARAMETERS
First, the ratio of actual to potentially-extractable soil water (SWR) is calculated for layer 1 (the upper layer) by the following equation:
where SW(1) is the initial soil water content of layer 1, LL(1) is the lower limit of plant-extractable soil water at layer 1, and DUL(1) = the drained upper limit of layer 1. SWR is then used to initialize the cumulative stage 2 soil evaporation (SUMES2).
If SWR is less than 0.9, then
and the cumulative stage 1 soil evaporation (SUMES1) is set equal to the upper limit of stage 1 soil evaporation (U). The time after the beginning of stage 2 evaporation (T) is then calculated.
If SWR is greater than 0.9, then both SUMES2 and T are set to 0 and the cumulative stage 1 soil evaporation (SUMES1) is calculated from the equation
These evaporation parameters come from the published model of Ritchie (1972).
The initiation subroutine then calculates the plant-extractable soil water (ESW(L)) for each layer
for use in output information and some calculations within the model.
The initialization program also calculates the cumulative depth of the soil profile (CUMDEP), the total soil water in the profile (TSW), the total plant-extractable soil water (TPESWP), the total soil water in the profile at the lower limit of plant extractable water (TLL), the total soil water in the profile at the drained upper limit (TDUL), and the total soil water in the profile at saturation (TSAT). These are used only to provide a printed output for users to quickly summarize the soil properties.
A weighting factor used to modify runoff curves numbers with near surface soil condition (WF(L)) is calculated from the equations
and
where XX is equal to 0 in the top layer, i.e. WF(L) = WX. In the other layers, XX is equal to the WX calculated for the layer above.
The intermediate variables used to calculate the run-off, CN1 and SMX, are calculated from the input curve number CN2.
Details of the runoff procedures are discussed in the next section.
INFILTRATION AND RUNOFF
Daily precipitation in mm is input into the model from the weather data file. If irrigation and/or precipitation occur on a day, the amount of irrigation (AIRR(J)) and precipitation (RAIN) are summed. The water balance subroutine calculates run-off by a modification of the USDA-Soil Conservation Service (SCS) curve number method (Williams et al. 1983). The SCS procedure uses the total precipitation from one or more storms which occur in a single day to estimate run-off, and excludes time as an explicit variable, i.e. rainfall intensity is ignored.
While the SCS procedure utilizes antecedent rainfall amounts to determine soil wetness and run-off, the procedure of Williams et al. for layered soils considers the wetness of the soil in the layers near the surface. This is accomplished by first calculating a sum (SUM), weighted for soil depth by the factor WF(L), of the relative amount of plant-extractable soil water in the profile. The SCS curve number retention parameter (R2) is then calculated from SUM and the maximum value of R2 (SMX):
where R2 cannot be less than 2.54 mm. The value SMX was calculated in the SOILRI subroutine.
Fig. 3.1 illustrates the SCS curve number concept with variations which allow for wet or dry conditions near the surface. The run-off curve concept is not expected to provide accurate run-off and infiltration information for a specific storm. The curve number concept was empirically derived to approximate the run-off volume when only daily rainfall is known. If a greater accuracy is required, a more physically based approach would be required. Such an approach would have to include information regarding storm intensities. Because storms vary both spatially and temporally, accurate modeling of infiltration and run-off would require more frequent rainfall measurements than daily values. These measurements would have to be taken in a dense network to provide spatially balanced input information. The knowledge of the right rainfall values for a particular site where the model is to be applied is one of the limiting factors affecting model accuracy.
Fig. 3.1. SCS curve number concept.
No run-off occurs if the temporary variable PB is less than zero, where
If run-off does occur, then infiltration is calculated as the difference between precipitation and run-off.
DRAINAGE
Water can be taken up by plants while drainage is occurring. Thus the drained upper limit soil water content is not always the appropriate upper limit of soil water availability. Many productive agricultural soils drain slowly, providing a possibly significant quantity of water to plants before drainage practically stops. In CERES-Wheat, the redistribution of water in the soil profile and drainage out of the rootzone are calculated using a functional model developed from field drainage information.
Using a cascading approach, water is moved downward from the top soil layer to lower layers. Drainage from a layer takes place when the soil water content (SW(L)) is between field saturation (SAT(L)) and the drained upper limit (DUL(L)).
For drainage calculations, the infiltration PINF is converted from mm to cm and a downward flux for each layer calculated (FLUX(L)). This information is needed for calculating leaching. When FLUX(L) is not equal to zero, the amount of water that the layer can hold (HOLD) between the current volumetric water content (SW(L)) and saturation (SAT(L)) is calculated.
If FLUX(L) is less than or equal to HOLD, an updated value of SW(L) is calculated prior to drainage.
If this new SW(L) is less than the drained upper limit of volumetric soil water in the layer (DUL(L)), no drainage occurs. If this SW(L) is greater than DUL(L), drainage (DRAIN) from the layer is calculated from SW(L), DUL(L), DLAYR(L), and SWCON, the whole profile drainage rate constant.
At this point, a drained value of SW(L) is calculated,
and a new value of FLUX(L), representing water moving into the layer below, is set equal to DRAIN.
If FLUX(L) is greater than HOLD, the water in excess of HOLD is passed directly to the layer below by saturated flow. The drainage is then calculated as follows:
An updated value for FLUX(L) is calculated
After calculating the water movement through all soil layers, any drainage from the bottom layer of the profile is equal to FLUX(L). For convenience, this value is converted to mm by multiplying by 10, and set equal to DRAIN. DRAIN then represents the total outflow from the lowest layer of the soil profile and is an available output variable for those interested in the time course of drainage out of the soil profile
EVAPOTRANSPIRATION
The soil water balance subroutine requires calculations for potential evaporation from the soil and plant surfaces. The equations to predict evaporation are primarily those used in the model of Ritchie (1972). The main difference between this part of the soil water balance subroutine and the Ritchie model is that a Priestly-Taylor (1972) equation for potential evapotranspiration is used instead of the Penman equation. This was done to eliminate the need for vapor pressure and wind inputs and under many circumstances it provides equal accuracy.
Calculation of potential evaporation requires an approximation of daytime temperature (TD) and the soil-plant reflection coefficient (ALBEDO) for solar radiation. For the approximation of the daytime temperature a weighted mean of the daily maximum (TEMPMX) and minimum (TEMPMN) air temperatures is used:
The combined crop and soil albedo (ALBEDO) is calculated from the model calculated leaf area index (LAI) and the input bare soil albedo (SALB). Prior to germination, ALBEDO is equal to SALB. For Stages 1 to 4 the value for ALBEDO is
For stages 4 to 6, ALBEDO is
An equilibrium evaporation rate (EEQ) defined in Priestly and Taylor (1972) is calculated from ALBEDO, TD, and the input solar radiation SOLRAD.
(Fig. 3.2). The units of EEQ is mm day-1 and SOLARD is MJ m-2 day-1. This empirical equation is a simplification of one containing long wave radiation approximations needed to calculate net radiation.
Fig. 3.2. EEQ/SOLARD versus TD.
The potential evaporation (EO) is calculated as the product of EEQ times 1.1. The constant 1.1 increases EEQ to a larger value to account for unsaturated air. This assumption is from the Priestley-Taylor assumption that the aerodynamic component of evaporation is proportional to the radiant energy into the system. The combustion equation assumes that the two terms are additive. When TEMPMX is greater than 24 C, the constant 1.1 is increased to account for advection
When TEMPMX is less than 5 C, the constant is reduced to account for cold temperature causing an additional decrease in EO due to stomata closure.
Fig. 3.3 graphically demonstrates how TEMPMX affects this constant multiplier for EEQ.
Fig. 3.3. TEMPMX relationships.
The potential rate of soil evaporation (EOS) is then calculated using the leaf area index, LAI. When LAI is less than 1,
and when LAI is greater than 1,
The calculation for the actual rate of soil evaporation (ES) is based on the assumption that there are two stages of soil evaporation. The first stage is limited by the energy available at the soil surface and continues until a soil-dependent upper limit is reached. This upper limit for stage 1 evaporation is expressed by the input U. After the upper limit of stage 1 is reached, soil evaporation enters stage 2. In stage 2, the rate of evaporation decreases proportionally to the time spent in stage 2.
The variables SUMES1 and SUMES2 are the sums of the soil evaporation (ES) in stages 1 and 2, respectively, and are used to determine which stage of soil evaporation is occurring during a day.
When rainfall or irrigation occurs during a day, and the infiltration into the upper layer is greater than or equal to SUMES1, SUMES1 is set back to 0. If, however, WINF is less than SUMES1, then SUMES1 is updated by the following equation:
Whenever SUMES1 is less than the upper limit of stage 1 evaporation (U), SUMES1 is updated daily by the following:
If the new value of SUMES1 is less than or equal to U, then
If, however, the new value of SUMES1 is greater than U, then
using the new value of SUMES1.
Should the value for SUMES1 exceed U, the soil evaporation enters stage 2, and SUMES2 is calculated.
The time after the beginning of stage 2 evaporation (T) is then calculated.
As the soil continues to dry during stage 2 evaporation, the value for T increases by 1 daily. The value for soil evaporation (ES) is calculated as follows:
When ES calculated in this manner is less the EOS, the value for ES is set to equal EOS.
When rainfall or irrigation occur during stage 2 evaporation and only slightly wet the soil, if WINF is less than SUMES2, then ES is set equal to the minimum of EOS, 0.8 * WINF, or ES + WINF.
During stage 2 evaporation, SUMES2 and T are updated daily.
If during Stage 2 evaporation, rainfall and/or irrigation wet the soil surface slightly more than described above such that WINF is greater than SUMES2, the value for T is reset to 0 and SUMES1 is calculated.
When this has occurred, the soil evaporation has re-entered Stage 1 and Stage 1 evaporation is calculated as described earlier.
After ES has been determined, the water is subtracted from the volumetric water content of soil layer 1 (SW(1)). If SW(1) decreases such that its value is less than SWEF * LL(1), ES is recalculated so that SW(1) does not become less than SWEF * LL(1). The term SWEF is considered the driest possible soil water content and is used to limit the drying of the soil by evaporation to a reasonable limit. The value of SWEF is approximately half the value of LL(1), with some variation depending on the layer thickness (DLAYR(1)).
The upward flow of water in the top 4 soil layers is next calculated to account for the movement of water, which usually moves toward the soil surface during evaporation. This upward flow is principally needed in the nitrogen subroutines to account for upward movement of nitrogen.
The variables THET1 and THET2 represent a normalized volumetric water content of layers L and L + 1, respectively as follows:
and values of THET1 and THET2 are constrained to be no less than 0. Thus the water content is normalized to the lower limit (LL(L) where the assumption is made that the diffusivity for all soils is similar. The assumed diffusivity (DBAR) is a function of the normalized water content for all soils
where DBAR is not allowed to be greater than 100. The unit of DBAR is cm day-1/2. The flow of water is then calculated:
A plot of DBAR, the soil water diffusivity, is shown in Fig. 3.4. The volumetric soil water in layers L and L + 1 are then increased and decreased by the amount of FLOW, and the new volumetric soil water is calculated.
Fig. 3.4. Soil water diffusivity.
For user information the cumulative soil evaporation after germination (CES) is updated daily,
The potential plant evaporation (EP) is calculated using simulated LAI values less than or equal to 3.
When LAI is greater than 3,
If EP + ES is greater than EO, then
Reducing EP from a potential value to an actual one requires the calculation of root water absorption.
ROOT WATER ABSORPTION
The root water absorption in CERES-Wheat is calculated using a law of the limiting approach whereby the soil resistance, the root resistance, or the atmospheric demand dominate the flow rate of water into the roots. Most of the details of this approach have been discussed by Ritchie (1984). The flow rates are calculated on the basis of water movement to a single root.
The maximum daily water uptake by roots in a layer (RWUMX) is assumed to be 0.03 cm3 of water per cm of root. This value sets the upper limit of water absorption by the roots as limited by axial root resistance. The potential root water uptake as influenced by soil water flow in a layer (RWU(L)) is
where RLV(L) is the root length density in the soil layer and ALOG is the FORTRAN equivalent of natural logarithm. This equation was derived from a radial flow to a single root and assumes that the hydraulic conductivity of all soils are similar when normalized to the lower limit value, LL(L). This assumption is more generally correct when the soil water content is near the lower limit. This equation also assumes that the water potential gradient between the root and the soil remains constant, even when the soil dries out. In reality, the water potential of the roots changes considerably throughout the day. However, because we are calculating daily values for water absorption, these less dynamic empiricisms provide sufficient detail for realistic uptake simulations.
If the calculated value for RWU(L) is greater than RWUMX, then RWU(L) is set equal to RWUMX (Fig. 3.5). The units of RWU(L) are then changed from cm3/cm of root to cm of water in a soil layer by the following:
and the total potential root water uptake from the entire root zone (TRWU) is calculated as the sum of RWU(L) for all soil layers with roots.
Fig. 3.5. The relationship used to calculate maximum root water absorption as related to O-' (the water content above the lower limit) and root length density (Lv). Also shown is the assumed maximum possible rate and the usual range of absorption when all the soil profile is at an optimum water content.
The terms TRWU and RWU(L) are potential water uptake values and represent the maximum possible water uptake from the soil profile and a root layer, respectively. A water uptake fraction (WUF) is calculated to reduce the potential root water uptake to the actual if necessary.
where EP1 is the potential transpiration in cm day-1. If WUF is less than or equal to 1.0, the plants are considered to be free of a water deficit so that the potential EP value is obtained. If WUF is greater than 1, then the actual plant evaporation is equal to the sum of the root absorption rates. Actual water uptake in each layer, RWU(L), is then calculated for the case when WUF is less than or equal to 1.0
Values of SW(L) are updated,
and the total soil water in the profile (TSW) is calculated:
The potentially extractable soil water in the entire soil profile (PESW) can now be calculated from TSW minus the total water content of the profile at the lower limit of the plant-extractable soil water (TLL):
The soil water deficits SWDF1 and SWDF2, used in the growth routine of CERES-Wheat, are now calculated:
SWDF1 and SWDF2 are constrained to be no greater than 1.0 (Fig. 3.6). If SWDF1 is not less than 1, the plant evaporation (EP) is assumed to be equal to the potential EP. If SWDF1 is less than 1, then the actual is equal to the maximum absorption rate (TRWU).
The 10. converts the units of EP to mm day-1. Actual daily value of soil plus plant evaporation now becomes
Fig. 3.6. Relationship used to calculate soil water factors, SWDF1 and SWDF2, to incorporate water stresses in the model.
To summarize the information on the soil water deficit influence on plants, two values are calculated: CSD1 and CSD2. These values are an average of the SWDF1 and SWDF2 during each growth stage. The values are not used in the model, but are given to provide information to the user for possible interpretation of the yield responses related to soil water deficits. The values appear on the printouts of the output summary sheet.
COLD HARDENING
Hardening, is a biochemical change that permits plants to survive under colder temperatures without being killed. Hardening occurs in the -0`C range (0 to -5`C or -10`C). These temperatures are not cold enough to kill the plants, but they are below the base temperature for them to be able to function. The longer a plant is in this state, the hardier it becomes and the longer it can withstand cold temperatures.
Plants go through two stages of hardening. The first stage is described above. The second stage occurs with even colder temperatures. A plant then can withstand the worst of conditions, though it can still be killed. (An unhardened plant will be killed by much higher temperatures than a hardened plant.)
The cold hardiness calculation in the CERES-Wheat model uses a hardiness index (HI) whose value varies between 0 and 2. Plants with HI values between 0 and 1 are considered to be in stage 1 hardening. Plants with HI values between 1 and 2 are in stage 2 hardening.
There are some limitations of this model for cold hardiness. It does not account for cold-sensitive wheat varieties planted where winter temperatures can kill the plants (we assumed only cold-resistant varieties would be planted in these areas). Genetic differences are not considered, as sufficient quantitative evidence to develop the necessary relationship for such differences is lacking. Also, this part of the model was not tested very thoroughly because we lacked the quantitative information from field experiments.
In the model, the vernalization process is started after the seed has germinated, because the activated seed can go through vernalization. In stage 9 and 1, the model automatically calls the cold subroutine. Hardening occurs from -1`C and -8`C. The quantitative expression in the model for the hardiness index is
where HI is the hardiness index and TEMPCR is the temperature at the crown. The HI is updated daily using this expression; HI in the left side of the equation comes from the previous day, and HI in the right side is the new HI. The expression is valid only between -1`C and 8`C. If the HI is 3.5, it would take 10 days for the HI to equal 1. Therefore, 10 days is the minimum time for a plant complete stage 1.
Stage 2 hardening is calculated when the first stage is greater than 1 and is then updated daily. (All plants must go through stage 1 before starting stage 2.) Any day the crown temperature is less than 0, the model takes the old value and adds 0.083 per day:
Under these conditions the plant will harden and HI will be at 2--maximum hardness--in 12 days.
Critical for the dehardening process is the maximum temperature (TEMPMX). Dehardening does not occur if TEMPMX is less than 10`C (Gusta and Flower 1976). If TEMPMX is 10, at 20`C, it loses 0.2 HI,
If the plant is in stage 2 hardening, the rate of daily dehardening is doubled.
After cold hardiness, leaf senescence is calculated if the minimum temperature is below -6`C. If these conditions exist, the model calculates the fraction of green leaf area (CR) killed using the expression
where SNOW is the snow depth in centimeters. The expression, empirically derived, accounts for slight modifications in senescence due to hardening and a strong influence related to the protection of the snow. Once CK is calculated the senescent leaf area (SENLA) is updated:
where PLA is the total plant leaf area developed during the plant's life.
If the plants have not been killed by the cold temperatures, a logical statement in the program prevents the green leaf area of each tiller from being less than 0.5 sq.cm. If the leaves and roots are senesced to that minimum point (0.5), the crown is assumed to contain 0.5 g of reserve carbohydrates for use in regrowth in much the same way. The stored carbohydrate is treated in the same way as the seed reserve was in helping to establish a new seedling before the plant becomes capable of independently producing enough assimilate for independent growth.
The next calculation is the death of tillers and plants. The threshold killing temperature (TEMKIL) for any death is calculated from the expression
This calculation shows that hardening has a lot to do with how much a plant will be killed by cold temperatures. TEMKIL is the threshold temperature by which plants begin to be killed. If a crown temperature is below that threshold, killing occurs by a reduction in tiller number if there is more than one tiller per plant. The expression for killing tillers is
If all but one tiller are killed, the plant population must be reduced. The reduction is done with an expression similar to tiller death:
If the plant population is reduced to less than five plants per square meter, the plant growth parts of the program are terminated and a message is printed stating that the crop has failed due to cold temperatures. To help users interpret possible reasons for yield reduction, a message is also written for days that the crop is damaged by cold temperatures.