Simulation of Wheat Growth and Development
This document is as a short documentation of the wheat crop-soil-model (CSM) ‘HumeWheat’ used in the current work of the agronomy group of CAU Kiel. The wheat specific part was developed within several projects with main contributions by Arne Ratjen (Ratjen et al., 2016, 2016; RATJEN and KAGE, 2015a; Ratjen and Kage, 2015), Ulf Böttcher, Adam Luig, Dorothee Neukam (Neukam et al., 2016) (Neukam et al., 2015) and others. It is in many aspects based on the CERES-Wheat model (CW) (RITCHIE, 2015), but uses different approaches especially for dry matter and nitrogen partitioning, leaf area expansion and dry matter production. Also the development part is extended in comparison to the original CERES wheat approach as it predicts BBCH stages on the basis of an explicit simulation of leaf initiation and appearance. ‘HumeWheat’ is a mechanistic CSM and runs on a daily time step. The model is of a modular design principle and was developed using an object oriented class component library (Kage and Stützel 1999a) implemented in the Borland Delphi Builder® V5 development environment. The model runs on daily time steps, but internal integration of the soil water balance module is on smaller time steps, the size dependent on the maximum changes of the soil water in the layers of the 1-D soil water balance model. Required external data inputs are the mean daily temperature [°C], precipitation [mm d-1], global radiation [W m-2], wind speed [m s-1], relative humidity [%] and values derived from these like vapor pressure [mbar] and saturation deficit [mbar].
1 Development
The phenological model module of HumeWheat is based on phenological Ceres Wheat (CW) model for winter wheat. It has been published originally in (Johnen et al., 2012a). The phenological process in CW is divided into nine developmental or growth stages (GS) according to (RITCHIE, 2015) (Table 1), the common approach in Europe, however is the BBCH stage (LANCASHIRE et al., 1991). Zadoks stages or equivalent BBCH values, however, are not calculated in the original phenology model of CERES-Wheat up to Version 3.0.
| GS | BBCH | Phase |
|---|---|---|
| 7 | 00 | Fallow |
| 8 | 00 | Sowing to germination |
| 9 | 05 | Germination to emergence |
| 1 | 10 | Emergence to terminal Spikelet initiation |
| 2 | 30 | Terminal spikelet to end of leaf growth and beginning of ear growth |
| 3 | 40 | End of leaf growth and beginning of ear growth to end of pre-anthesis ear growth |
| 4 | 57 | End of pre-anthesis ear growth to beginning of grain filling |
| 5 | 71 | Grain filling |
| 6 | 90 | End of grain filling to harvest |
The CW model mainly describes phenological development through the change of GS over time. In CWm we added an explicit simulation of three additional state variables, the leaf number on the main stem (nLMS), the number of initiated leaves on the main stem (inLMS) and the BBCH stage. At emergence nLMS is set to one and inLMS is set to 4 (Brooking and Jamieson, 2002). This is because the embryo at the time of sowing has already three leaf primordia (KIRBY et al., 1985) and that an additional primordium is assumed to be produced until emergence. The rate of change of number of initiated leaves is calculated according to Eq. 1 from emergence until the last leaf primordium is formed:
\frac{d\,in{{L}_{MS}}}{dt}=\frac{\max (0,T-{{T}_{b}})}{Plast} \tag{1}
Where Plast is the plastochron (°Cd) T is the air temperature and Tb is the base temperature. The formation of leaf primordia which produce leaves is finished some time before the double ridge stage (Bonnet, 1935; Bonnet, 1936; Muschik, 1957). The time of the cessation of leaf primordia initiation, however, is not exactly coupled to a certain GS in the CW model and is also not linked to a certain BBCH stage. From available data of BBCH 37, i.e. the time when the last leaf appears, it was possible to estimate an approximate GS corresponding to the formation of the final leaf primordium (GS_{flp}) ex post. The number of visible leaves on the main stem nLMS is calculated accordingly (Eq. 2):
\frac{d\,n{{L}_{MS}}}{dt}=\frac{\max (0,T-{{T}_{b}})}{Phyll} \tag{2}
where Phyll denotes the phyllochron (°Cd).
1.1 Calculation of changes of GS and BBCH
The main state variables of the phenological model are the GS-stage according to Ritchie and the BBCH-stage more commonly used in the description of wheat phenology. Their rate of change is calculated stage specific.
1.1.1 Sowing to emergence (GS 8 and 9)
We simplified the model by pooling together the GS stages eight and nine, which represent the germination and emergence processes, because no data were available to distinguish between both events. The redefined parameter P9, which now governs the rate of development during the time from sowing to emergence, was calculated from observed thermal times from sowing to emergence .
\frac{dGS}{dt}=\frac{\max (0,T-{{T}_{b}})}{P9} \tag{3}
The rate of change for BBCH is ten time higher as for GS as the event of crop emergence is equivalent to 10.
\frac{dBBCH}{dt}=10\frac{\max (0,T-{{T}_{b}})}{P9} \tag{4}
1.1.2 Emergence to terminal Spikelet initiation (GS 1)
Temperature, vernalization, photoperiod, and phyllochron interval determine the development rate during GS 1. According to CW the lowest value of two functions, depending on the vernalization stage and the actual photoperiod in interaction with the actual effective temperature limits the development rate (Eq. 5).
\frac{dGS}{dt}=\frac{\max (0,T-{{T}_{b}})\cdot \min (f(V),f(P))}{(400\cdot Phyll\ /\ 95)} \tag{5}
The functions f(V) and f(P) as well as the Phyllocron (Phyll) are genotype specific.
1.1.2.1 Vernalisation
The effect of vernalisation on the development rate is calculated from the accumulated sum of “vernalisation days”. The rate of change of the vernalisation days is calculated from the daily mean temperature and a vernalisation response function which is characterized by 4 cardinal temperatures (Figure 1).
The vernalisation rate is acccumulated in the state variable VernalisationDays.
\frac{dV}{dt}=\max (0, fV) \tag{6}
The value of the f(V) in Eq. 5 is calculated from the actual value of accumulated number of vernalisation dates according to:
f(V) = 1-k\cdot \left(50-V \right) \tag{7}
The values of k thereby is genotype specific and according to Ritchie et al. (1998) is coded in the parameter value P1V. The value of 50 in Eq. 7 is the maximum number of vernalisation days. The vernalisation days are reset to zero after the maximum number of vernalisation days is reached.
The value k is the slope of the linear relation between the vernalisation days and the development rate. It can be calculated from the parameter P1V according to
k = \frac {P1V+0.55}{183} \tag{8} For a larger number of wheat cultivars in Germany (Johnen et al., 2012b) found a value of P1V of 2.84. Higher values of P1V indicate a higher vernalisation requirement of the cultivar, i.e. a lower development rate at low vernalisation days which is typical for cultivars of continental origin. A low value of P1V indicates a low vernalisation requirement, typical for summer crops of wheat or for wheat cultivars grown in Mediterranean regions.
1.1.2.2 Photoperiod
The effect of photoperiod on the development rate (f(P)) is calculated from the actual photoperiod and a photoperiod response function
f(P)=1-C \cdot\left(20-P\right)^2 \tag{9}
The parameter C is coded in the parameter P1D. The value of 20 in Eq. 9 is the maximum photoperiod. The value of C is the slope of the quadratic relation between the photoperiod and the development rate. It can be calculated from the parameter P1D according to
C=P1D/500 \tag{10}
Typical values of P1D are in the range of 1 to 3. Higher values of P1D indicate a higher photoperiod requirement of the cultivar, i.e. a lower development rate at low photoperiods which is typical for cultivars of northern origin. A low value of P1D indicates a low photoperiod requirement, typical for summer crops of wheat or for wheat cultivars grown in mediterranean regions. In the study of (Johnen et al., 2012c) a value of P1D of 2.76 was found as a mean for a larger number of wheat cultivars in Germany.
1.1.2.3 Leaf number and leaf initiation
At early development stages (BBCH<30) BBCH stages are determined by the number of leaves and tillers present. Their rate of change is determined by the Phyllochron (Phyll), i.e. the inverse of the phyllochron is the emergence rate of leaves and tillers per effective day temperature (@eq-dBBCH_dt_GS1).
\frac{dBBCH}{dt}=\frac{\max (0,T-{{T}_{b}})}{Phyll} \tag{11}
At BBCH 13.5 there is a switch to BBCH 21, because simultaneously to the appearance of the fourth leaf the first tiller is emerging and mainstem leaf appearance is then associated with tiller appearance until BBCH 30.
The plastochron is the inverse of the rate of leaf initiation. It is assumed to be constant during the time from sowing to the end of leaf initation. It defines the rate of leaf initiation and together with the initial number of leaf primordia of the embryo at sowing the number of initiated leaves and with the GS stage when leaf initation ceases the number of leaves on the main stem.
If \ \ BBCH >= 13.5 \ and \ BBCH<20; \ BBCH = BBCH+7.5 \tag{12}
1.1.3 Terminal spikelet to end of leaf growth and beginning of ear growth (GS 2)
Within the model the implicit and simplifying assumption is made that GS 2 stage is more or less closely associated with BBCH 30 (Baker and Gallagher, 1983), but that the initiation of the terminal leaf primordium may be reached a bit earlier at a stage we named GS_{flp}.
The GS change rate between GS 2 and 3 is calculated as the ratio of the daily temperature minus the base temperature (Tb, assumed to be 0°C) to the thermal time interval GS 2 to GS 3. This thermal time interval is calculated by the number of leaves not yet expressed at GS 2 minus two (fL), multiplied by the phyllochron plus the temperature sum between BBCH 37 and 39, i.e. ‘PH39’:.
\frac{dGS}{dt}=\frac{\max (0,T-{{T}_{b}})}{fL\cdot Phyll+Ph39} \tag{13}
thereby the term fL denotes the leaves which have to appear after GS 2 is reached, calculated from the number of initiated leaves at GS 2 minus 2. The substraction of two leaf primordia was made in accordance to the assumption that at least the last two vegetative primordia may not produce visible leaves, they are therefore sometimes called labile primordia (GRIFFITHS et al., 1985).
In the CWm model we therefore calculated a rate of node emission and introduced a parameter we called ‘TSumInternode’ analogous to the phyllochron. The numerical value of this parameter is the inverse of the slope of the relation of node numbers to the temperature sum from BBCH 30. The rate of change of BBCH stages is calculated until BBCH 37 (final leaf begins to emerge) from the ratio of the effective day temperature and TSumInternode (Eq. 11).
The BBCH stage 37 is reached when the number of the emerged leaves is equal to the number of leaf primordia minus two. The subtraction of two leaf primordia was made in accordance to the assumption that at least the last two vegetative primordia may not produce visible leaves, they are therefore sometimes called labile primordia (GRIFFITHS et al., 1985).
The length of the thermal time from sowing to the initiation of the last leaf primordium together with the plastochron now determines the number of initiated leaves. This number is variable and the new algorithm of CWm therefore also predicts a variable length of GS 2 in terms of temperature sum. In CWm the rate of development between BBCH 37 and BBCH 39, however, is again described by a certain temperature sum expressed as the newly added parameter ‘PH39’ (Eq. 15). It was introduced because according to our data analysis the thermal time between BBCH 37 and 39 ratings was significantly longer than one phyllochron.
\frac{dBBCH}{dt}=\frac{\max (0,T-{{T}_{b}})}{\text{TSumInternode}} \tag{14}
If nL_MS > inL_MS-2 and If BBCH< 37 then BBCH = 37 If BBCH>=37:
\frac{dBBCH}{dt}= \min \left( \frac{2\max (0,T-{{T}_{b}})}{\text{Ph39}},(40-BBCH) \right) \tag{15}
1.1.4 End of leaf growth and beginning of ear growth to end of pre-anthesis ear growth (GS3)
The period from the end of leaf growth to the beginning of ear growth is for the GS development rate defined by a length of 2 phytochrons
\frac{dGS}{dt}=\frac{\max (0,T-{{T}_{b}})}{2Phyll} \tag{16}
The rate of change of BBCH stages is during that stage simply synchronized with the GS rate according to equation Eq. 17.
\frac{dBBCH}{dt}=(4+1.7(GS-3))\cdot 10-BBCH \tag{17}
1.1.5 End of pre-anthesis ear growth to beginning of grain filling (GS4)
The period from the end of pre-anthesis ear growth to the beginning of grain filling is for the GS development rate defined by a length of 200 degree days:
\frac{dGS}{dt}=\frac{\max (0,T-{{T}_{b}})}{200} \tag{18}
The rate of change of BBCH stages is during that stage simply synchronized with the GS rate according to equation Eq. 19.
\frac{dBBCH}{dt}=(5.7+1.4(GS-4))\cdot 10-BBCH \tag{19}
1.1.6 Grain filling (GS5)
The period of grain filling is for the GS development rate defined by a genotype specific temperature sum and the effective temperature during grain filling. The latter is calculated as the mean temperature minus the base temperature. The effective temperature is divided by the temperature sum for grain filling (TSumGF) to get the rate of change of GS 5:
\frac{dGS}{dt}=\frac{\max (0,(T-{{T}_{b}})-1)}{TSumGF} \tag{20}
TSum_{GF} = \frac{(P5+21.5)}{0.05} \tag{21}
The rate of change of BBCH stages is during that stage simply synchronized with the GS rate according to equation Eq. 22.
\frac{dBBCH}{dt}=(7.1+1.9(GS-5))\cdot 10-BBCH \tag{22}
1.1.7 End of grain filling to harvest (GS6)
The period from the end of grain filling to harvest is for the GS development rate defined by a length of 250 degree days:
\frac{dGS}{dt}=\frac{\max (0,T-{{T}_{b}})}{250} \tag{23}
The rate of change of BBCH stages is during that stage simply synchronized with the GS rate according to equation Eq. 24.
\frac{dBBCH}{dt}=(9+(GS-6))\cdot 10-BBCH \tag{24}
2 Dry matter production
The daily dry matter accumulation (\frac{d{W_{tot}}}{dt}) of the canopy is calculated according to the LUE approach, using the simple big-leaf canopy model. Within the source code the term carbo is used to denote the daily growth rate. The potential daily dry matter accumulation rate is a function of intercepted PAR multiplied with the potential light use efficiency LUEpot). LUE is affected by extremes of daily mean temperature, water limitation, nitrogen concentration of the canopy and the external CO2 concentration of the atmosphere as described in the following sections.
flowchart LR
A(LUE<sub>eff</sup>) -.-> F
C(f<sub>T</sup>) -.-> A
D(SWDF) -.-> A
E(Q) --> F{dW/dt}
G(I) --> E
H(f<sub>int</sub>) -.-> E
I(LAI) -.-> H
J((LUE<sub>pot</sub>)) -.-> A
K(f<sub>SLN</sub>) -.-> A
L(SLN) -.-> K
M(f<sub>CO2</sub>) -.-> A
N(CO<sub>2</sub>) -.-> M
O(CWSI) -.-> M
2.1 Radiation interception
According to the Lambert-Beer law the fraction of the intercepted PAR is a function of leaf area index (LAI) and extinction coefficient (kPAR). Canopy reflection within the PAR range is thereby lumped into the empirical kPAR value. The state variable LAI is calculated in the leaf area module. The fraction of intercepted PAR (fint) is calculated as follows:
Q= PAR\text{(t)} \cdot \left ( 1-{e^{-{k_{PAR}} \cdot LAI\text{(t)}}} \right) \tag{25} This is a very simplified approach, not considering fractions of direct and diffuse radiation and their complex interactions with the canopy. To account for the effect of LAI on the extinction coefficient, the extinction coefficient is calculated as a function of LAI and the maximum extinction coefficient (kPARini) according to the following equation:
k_{PAR} = min \left( k_{PAR_{max}},\ max(k_{ini}+ln(LAI) \cdot k_{dec}, \ k_{PAR}) \right) \tag{26}
where kini is the initial extinction coefficient and kdec is the extinction coefficient decrease per unit of LAI. The standard extinction coefficient (kPARmax) is set to 0.7 in the model.
2.2 Limiting factors on LUE
Suboptimal temperature and drought stress reduces dry matter production rates below the potential. The influence of suboptimal conditions is quantified by a temperature (fT) and a soil water deficit factor (SWDF).
2.2.1 Temperature
Daily mean temperature is calculated in a classical way from observed maximum and minimum temperature. if difference between min and max temperature is high, the maximum temperature is weighted stronger.
if \left((T_{max} - T_{min}) > \Delta T_{crit}\right) then \\ \quad T_{mean} = T_{max} \cdot T_{max_{wF}} + T_{min} \cdot (1 - T_{max_{wF}}) \tag{27}
The weighting factor T_{max_{wF}} is taken as 0.567. The critical temperature difference \Delta T_{crit} is set to 8°C.
The effect of temperature on LUE is calculated using a temperature effect factor fT [0..1] from the daily mean temperature. The base for its calculation is the equation of (Wang and Engel, 2002) (Eq. 28, Figure 6).
fT = \frac{ \left( 2 \cdot (T - T_{min})^{\alpha} \cdot (T_{opt} - T_{min})^{\alpha} - (T - T_{min})^{(2 \cdot \alpha)} \right)} { (T_{opt} - T_{min})^{(2 \cdot \alpha)}} \tag{28}
with \alpha = \frac{ ln(2)}{ln((T_{max}-T_{min})/(T_{opt}-T_{min}))} \tag{29} For Figure 6 the parameters Tmin = 0°C, Topt = 20°C and Tmax = 35°C were used.
If available an canopy surface temperature can be used to calculate the temperature effect factor. The minimum of the obtained fT value is used for the calculation of the daily dry matter accumulation rate.
2.2.2 Drought stress
For stomatal aperture, the diffusion paths and therefore the resistances/conductivities for CO2 and water evaporation are not identical, because of the liquid-phase components of the CO2 pathway into the mesophyll. Consequently, a non-linear relation between transpiration and LUE reduction is assumed, according to (Ferreyra et al., 2003). The formulation based on the RT factor technically calculated in the soil water module: RT is defined as the ratio between the sum of Tact and interception evaporation (Ip) divided by the sum of Tpot and Ip. The factor is calculated as the ratio
RT= \frac{ (T_{act}+I_p)}{(T_{pot}+I_p)} \tag{30}
Within the source code of HUMEWheat the factor RT is called TransIntRatio. The soil water deficit factor (SWDF) is calculated from the RT as follows:
SWDF = 1-( 1-RT )^{pSWDF} \tag{31}
The value of pSWDF thereby defines the non-linearity of the relation between RT and SWDF. The value of pSWDF is set to 2.5 in the model. Figure 7 shows the relation between RT and SWDF for different values of pSWDF.
2.2.3 Nitrogen limitations
Consistent with Sinclair and Horie (1989), Meinke et al. (1997) and Meinke et al. (1998a) HumeWheat uses an SLN threshold value SLNcrit [g N.cm-2] below which a LUE reduction factor (fSLN) is applied.
In a first step the optimum SLN value (SLNcrit) is calculated from the initial SLN value (SLN_{crit_{int}}) and the green area index (GAI) according to the following equation:
SLN_{crit} =SLN_{crit_{int}}+SLN_{crit_{inc}} \cdot GAI; \tag{32}
The values of the parameters SLN_{crit_{int}} and SLN_{crit_{inc}} are set to 3.74 and -0.228, respectively (RATJEN and KAGE, 2015b). The value of SLN_{crit_{int}} is the initial SLN value at GAI = 0. The value of SLN_{crit_{inc}} is the increase of the SLNcrit value per unit of GAI.
The Figure 9 from RATJEN and KAGE (2015b) gives the original source of this relationship derived from an analysis of experimental data using quartile regression.
Before booting SLNcrit is always lower than SLNopt.
A nitrogen limitation factor SLNI is calculated from the ratio of the actual SLN value (SLN) and the critical SLN value (SLNcrit) according to the following equation:
SLNI = min \left( 1, \frac {SLN} {SLN_{crit}} \right) \tag{33}
The obtained value then SLNI is used to calculate the nitrogen limitation factor Carbored according to the following non-linear equation:
Carbo_{red} = min\left(1, max \left(0, SLNI_a. + SLNI_b \cdot SLNI + SLNI_c \cdot SLNI^2 \right)\right) \tag{34}
The values of the parameters SLNIa, SLNIb and SLNIc are set to -0.197, 2.8 and -1.6, respectively (RATJEN and KAGE, 2015b). The value of SLNIa is the intercept of the Carbored function. The value of SLNIb is the slope of the Carbored function. The value of SLNIc is the curvature of the Carbored function. The Carbored function is shown in Figure 11.
The name of the variable Carbored is derived from the original CERES-Wheat source code where it is primarily used to calculate the reduction of the photosynthesis during ripening/senescence. We use it here to calculate the general reduction of the LUE due to nitrogen limitation.
2.2.4 CO2 concentration
The effect of a varying CO2 concentration on LUE is optionally calculated using a CO2 effect factor fCO2~ [0..1]. Thereby an interaction with the crop water stress status is assumed. The CO2 effect factor is calculated as follows:
f_{CO2} = max\left( f_{CO_2min}, \ f_{CO_2min} \cdot \left(\left[CO_2\right]-\tau \right)^{fCO_2} \right) \tag{35}
with
f_{CO_2min} = (f_{CO_2}scale+CWSI \cdot f_{CWSI}) \tag{36}
The crop water stress index is calculated from 1-TransIntRatio (0 = no drought stress, 1 = maximum drought stress). The interaction of CWSI and CO2 concentration is described by the parameter fCWSI. The parameter fCO2~scale is the CO2 effect factor at the CO2 concentration of the compensation point (\tau). The parameter fCO2~ is the slope of the CO2 effect factor curve. The parameter \tau is the CO2 concentration at which the CO2 effect factor is equal to fCO2~scale. The CO2 effect factor is limited to a value between 0 and 1. It follows from Figure 12 that the CO2 effect factor increases with increasing CO2 concentration and increasing CWSI, which means that the effect of CO2 on LUE is stronger under drought stress conditions.
2.3 Daily dry matter accumulation
The daily dry matter accumulation rate of the canopy is calculated as follows:
PCARB = LUE_{pot} \cdot PAR \cdot f_{int} \cdot min(f_T,Tempf_{surface}) \cdot f_{CO2} \cdot min(SWDF, Carbo_{red}) \tag{37}
where LUE is the light use efficiency, PAR is the photosynthetically active radiation, f_{int} is the fraction of intercepted PAR, f_T is the temperature effect factor, Tempf_{surface} is the temperature effect factor calculated from the canopy surface temperature, f_{CO2} is the CO2 effect factor, SWDF is the soil water deficit factor and Carbo_{red} is the nitrogen limitation factor.
The daily dry matter per plant is calculated as follows:
CARBO = PCARB / Plants \tag{38}
From the daily dry matter accumulation rate and the amount of intercepted radiation and the effective radiation use efficiency the daily dry matter accumulation rate of the canopy is calculated.
LUE_{eff} = \frac {PCARB}{IPAR} \tag{39}
where IPAR is the intercepted PAR. The effective radiation use efficiency thereby sums up the different limitations of the LUE into a single value.
3 Dry Matter Partitioning
The daily dry matter accumulation rate of the canopy is partitioned into the different plant organs. The partitioning of the daily dry matter accumulation rate into the different plant organs is in contrast to many other crop growth models based on the concept of allometric relationships between the organs of the crop. The other governing principle is the concept N dilution, i.e. decreasing N concentrations of the organs.
Effects of nitrogen shortage and drought stress on the partitioning of dry matter are considered in the model, but can be switched off by setting the respectives model options.
3.1 Initialisation
The initial leaf weight per plant after emergence is calculated from the initial leaf area and the specific leaf area (SLA) according to the following equation:
LFWT_{pl} = LA_{ini} / SLA_{pot} \tag{40}
The stem weight per plant then is calculated from an allometric relationship between leaf and stem weight according to the following equation:
STMWT_{pl} = exp(g * ln(LFWT_{pl} * plants) + h) / plants \tag{41}
where g and h are parameters of the allometric relationship. The default values of g and h are set to 1.46 and -2.13 in the model.
The seed reserves from the parameter SEEDRVini:
SEEDRV = SEEDRV_{ini} \tag{42}
The nitrogen amount in the leaf fraction is calculated from the leaf nitrogen concentration under optimum N supply and the initial leaf weight per plant according to the following equation:
NLeaf_{pl} = LFWT_{pl} * Nc_{opt_{Leaf}} / 100 \tag{43}
Accordingly, the nitrogen amount in the stem fraction is calculated from the stem nitrogen concentration under optimum N supply and the initial stem weight per plant according to the following equation:
NStem_{pl} = STMWT_{pl} * NoptStem / 100 \tag{44}
The initial amount of N in the shoot of a single wheat plant is calculated as the sum of the N in the leaf and stem fraction:
NShoot_{pl} = NStoragepool_{pl} + NLeaf_{pl} + NStem_{pl} + NGrain_{pl} \tag{45}
The N amount in the root fraction is calculated from the initial shoot N amount and a fraction of the initial shoot N amount according to the following equation : NRoot_{pl} = NShoot_{pl} * fRootN_{ini} \tag{46}
3.2 SeedRV change
The seed reserves contribute to the daily dry matter accumulation rate of the canopy. Their rate of change is calculated as follows:
\frac{SEEDRV}{dt} = -k_{SEEDRV} * SEEDRV \tag{47}
where k_{SEEDRV} is the rate of change of seed reserves. The value of k_{SEEDRV} is set to 0.15 in the model, i.e. the daily change rate is 15% of the actual amount.
The rate of change of seed reserves is added to the assimilate flow calles ‘Assiflow’ in the model:
$$Assiflow = Assiflow - SEEDRV.c;
3.3 Allometric Leaf Stem Partitioning
3.3.1 Drought effects
The partitioning of dry matter between leaf and stem is influenced by the soil water potential. Under water stress conditions the partitioning of dry matter to the stem is increased (Figure 13).
Within the dry matter module the variable kfd describes the influence of psiroot, i.e. the absolute, non negative weighted average of the soil hydraulic potential on leaf stem partitioning.
The value of k_f is calculated as follows:
kf_d = max(1,\ 1 + (| psi_{root} | - psi_{crit}) \cdot psi_s) \tag{48}
where psiroot is the soil hydraulic potential, psicrit is the critical soil water potential at which the leaf stem partitioning is 1 and psis is the slope of the relationship between leaf stem partitioning and soil water potential. The default value of psis is set to 0.296 in the model. The value of psicrit is set to a default of 2.24 MPa in the model. These parameters are from Ratjen et al. (2016).
The value of psicrit is the critical soil water potential at which the leaf stem partitioning is 1. The value of psis is the slope of the relationship between leaf stem partitioning and soil water potential. The default value of psis is set to 0.296 in the model. The value of psicrit is set to a default of 2.24 MPa in the model.
3.3.2 Nitrogen effects
Under nitrogen deficiency a higher fraction of dry matter is partitioned to the stem. The partitioning of dry matter between leaf and stem is calculated from the nitrogen nutrition index (Eq. 49).
kf_n = max(1,\ 1 + (NNI_{crit} - NNI) \cdot NNI_{inc}) \tag{49}
where NNIcrit is the critical nitrogen nutrition index at which the leaf stem partitioning is 1, NNIinc is the increase of the leaf stem partitioning per unit of NNI and NNI is the actual nitrogen nutrition index. The default value of NNIinc is set to 0.595 in the model. The value of NNIcrit is set to a default of 1.377 in the model.
The stem fraction of shoot dry matter increase is calculated from an allometric relation ship:
f_{stem*} = 1 - \frac{1}{ 1 + e^h \cdot LFWT^{g-1} \cdot g} \tag{50}
where LFWT is the leaf DM per square meter and h and g are parameters of the allometric relationship. The default values of g and h are set to 1.46 and -2.13 in the model.
f_{stem} = min(1, \ f_{Stem*} \cdot max(kf_d, \ kf_n))
3.4 Fine root growth
The fraction of assimilates allocated to fine roots is calculated from a linear decrease of fine root mass with increasing temperature sum.
f_{fineroot} = max \left( 0, \left(f_{FineRoot0} - f_{FineRootDec} \cdot Tempsum \right) \right) \tag{51}
where f_{FineRoot0} is the initial fraction of assimilates allocated to fine roots, f_{FineRootDec} is the decrease of the fraction of assimilates allocated to fine roots per unit of temperature sum and Tempsum is the temperature sum since sowing. The default values of f_{FineRoot0} and f_{FineRootDec} are set to 0.653 and 0.000501 in the model.
Historically the fraction of assimilated allocated to the shoot is called PTF in Ceres Wheat.
PTF = min \left( 1, 1 - f_{fineroot} \right)
Effects of drought stress and nitrogen shortage on the partitioning of dry matter between shoot and root are considered in the model, but can be switched off by setting the respectives model options. These effects are implemented by altering the value of PTF and ffineroot is re-calculated from PTF.
Under drought stress the initial value of PTF is decreased by 0.1 times the difference between 1 and the ratio of actual to potential tranpiration (TransIntRatio) according to the following equation:
PTF = PTF - 0.1 \cdot (1 - TransIntRatio) \tag{52}
This is a preliminary parameterization of the effect of drought stress on the partitioning of dry matter between shoot and root. The value of 0.1 is a first guess and should be further calibrated.
Under nitrogen shortage the initial value of PTF is decreased by 0.1 times the difference between 1 and the nitrogen nutrition index (NNI) according to the following equation:
PTF = PTF - 0.1 \cdot (1 - NNI) \tag{53}
This is a preliminary parameterization of the effect of nitrogen shortage on the partitioning of dry matter between shoot and root. The value of 0.1 is a first guess and should be further calibrated.
The fraction of assimilates allocated to fine roots is limited to a value between 0 and 1 and is used to calculate the daily change of the fine root mass per plant:
f_{fineroot} = max(0, 1 - PTF) \tag{54}
3.5 Leaf nitrogen concentration
After the BBCH stage of the crop reaches the value EC_{crit_{NcLeaf}}, which is set to 27, i.e. before shooting, the leaf nitrogen concentration [% N of DM] is calculated as follows:
Nc_{opt_{leaf}} = Nc_{leaf_{Vf1}} \cdot LFWT + Nc_{leaf_{Vf2}} \tag{55} thereby Nc_{leaf_{Vf1}} and Nc_{leaf_{Vf2}} are parameters of the linear relationship between leaf nitrogen concentration and leaf dry matter. The default values of Nc_{leaf_{Vf1}} and Nc_{leaf_{Vf2}} are set to -0.0061 and 5.9543 in the model.
In a transitional phase between autumn and spring, the leaf nitrogen concentration is calculated as follows:
\begin{align*} \frac{ Nc_{Leaf_{Winter}}}{dt} & = (Nc_{leaf_{winter}} - Nc_{leaf_{min}}) \cdot rgr_{NcLeafWinter} \cdot \\ & GlobRad \cdot \frac{ 1 - (Nc_{leaf_{winter}}- Nc_{leaf_{min}})} {Nc_{Leaf_{Vf1}} \cdot LFWT + Nc_{Leaf_{Vf2}} - Nc_{leaf_{min}}} \end{align*} \tag{56}
where Nc_{leaf_{min}} is the minimum leaf nitrogen concentration, rgr_{NcLeafWinter} is the relative growth rate of the leaf nitrogen concentration, and GlobRad is the global radiation. The default values of Nc_{leaf_{min}} and rgr_{NcLeafWinter} are set to 1.5 and 0.0001 in the model. The value of Nc_{leaf_{min}} is the minimum leaf nitrogen concentration. The value of rgr_{NcLeafWinter} is the relative growth rate of the leaf nitrogen concentration. The value of GlobRad is the global radiation. The value of Nc_{leaf_{min}} is the minimum leaf nitrogen concentration. The value of rgr_{NcLeafWinter} is the relative growth rate of the leaf nitrogen concentration. The value of GlobRad is the global radiation.
if DayOfYear < 150 then if TEMPM > 0 then // transitional phase between autumn and spring
Nc_optLeaf = NcLeafWinter \tag{57}
else
3.5.1 Calc_Leaf_N_Change
Above the critical BBCH stage (EC_{crit_{NcLeaf}}, set to EC 27, i.e. end of tillering) the optimum leaf nitrogen concentration is calculated from a linear dilution curve:
Nc_{opt_{Leaf}} = Nc_{Leaf_{Vf1}} \cdot LFWT_{m2} + Nc_{Leaf_{Vf2}}
where Nc_{Leaf_{Vf1}} and Nc_{Leaf_{Vf2}} are parameters of the linear relationship between leaf nitrogen concentration and leaf dry matter. The default values of Nc_{Leaf_{Vf1}} and Nc_{Leaf_{Vf2}} are set to -0.0061 and 5.9543 in the model.
The daily change of the leaf nitrogen amount per plant is calculated from the daily change of leaf dry matter and the leaf nitrogen concentration under optimal N supply and the difference between the actual leaf nitrogen amount and the optimal leaf nitrogen amount in the leaf fraction of a plant:
\frac {N_{Leaf_{pl}}}{dt} = \frac {LFWT_{pl}}{dt} \cdot Nc_{opt_{Leaf}} / 100 + Nc_{opt_{Leaf}} / 100 \cdot LFWT_{pl} - N_{Leaf_{pl}} \tag{58}
3.5.2 Stem N
A similar approach is used to calculate the stem nitrogen concentration under optimal N supply. The stem nitrogen concentration under optimal N supply is, however, calculated from a non-linear dilution curve:
NearOptNCStem = 0.98 \cdot (1 / (NcStem_a + NcStem_b \cdot STMWT_{m2}))
NcStem <- function(STMWT, NcStem_a, NcStem_b){
NcStem <- 0.98 * (1 / (NcStem_a + NcStem_b * STMWT))
return(NcStem)
}
STMWT <- seq(0, 400, 1)
NcStem_a <- 0.1475
NcStem_b <- 0.0011
df_NcStem <- data.frame(STMWT = STMWT)
df_NcStem$NcStem <- NcStem(df_NcStem$STMWT, NcStem_a, NcStem_b)
df_NcStem$NcStemNearOpt <- 0.98*df_NcStem$NcStem
ggplot(data = df_NcStem) +
geom_line(aes(x = STMWT, y = NcStem, color=NcStem),size=2) +
geom_line(aes(x = STMWT, y = NcStemNearOpt), size=1, linetype=2) +
labs(x = "Stem dry matter [g m^-2]", y = "NcStem [%]") +
theme_bw(base_size = 20) +
scale_color_gradient(low = "yellow", high = "darkgreen") +
scale_y_continuous(limits = c(0, 5)) 
3.6 Grain filling
3.6.1 Grain number
The model for grain number follows the assumption that the grain number per square meter (GPSM) is a function of the shoot dry matter at mid flowering i.e. EC65 (DM~65), the nitrogen nutrition index at EC 60 (NNI~60), and the average photothermal quotient from 45 d preceding anthesis (Q45):
GPSM= G4\cdot ln \left( DM_{65} \cdot min \left( NNI_{60}\right) \cdot Q45 \right)^{a} \tag{59}
where GM4 is a parameter of the regression equation. The default values of G4 and a are set to149.153 and 2.6233 in the model (see Figure 17). Within the source code the parameter G4 ist called GM4 and the parameter a is named GM4_a. See (Ratjen et al., 2012) for the derivation of this approach .
For the estimation of grains per plant in order to calculate the initial value of grain N an early estimate of GPSM is derived from the shoot dry matter at anthesis (TOPWT_m2) and the photothermal quotient 45 d preceding anthesis (Q45):
if (EC >= 40) and (NNI60 <= 0) and if ln(TOPWT_{m2} \cdot Q45) > 0 then
GPPVAR = \left ( GM4 \cdot \left ( ln \left ( TOPWT_{m2} \cdot Q45 \right) \right)^{GM4_2} \right) / plants \tag{60}
where GM4~2 is a parameter of the regression equation. The default value of GM4~2 is set to 0.0001 in the model.
3.6.2 Initialisation
The initial change of the grain weight is calculated using a sliding N and C accumulation of initial grain N & C from EC 40 to 65 from the parameter GRNWT_{ini} and the missing fraction of that weight. Initial grain weight per plant and the grain filling rate according to the following equation:
\frac{GRNWT_{pl}}{dt} = GRNWT_{ini} / 1000 \cdot GPPVAR \cdot min\left( 1, \frac{(EC - 40)}{ (65 - 40)}\right) - GRNWT_{pl} \tag{61}
where GRNWT_{ini} is the initial grain weight per plant, GPPVAR is the grain filling rate, EC is the EC stage. The value of GRNWT_{ini} is set to 0.1 in the model. The value of GPPVAR is set to 0.0001 in the model.
Because the initial weight comes from the stem fraction, the rate of change of DM stem has to be adjusted:
\frac{STMWT_{pl}}{dt} = \frac{STMWT_{pl}}{dt} - \frac{GRNWT_{pl}}{dt} \tag{62}
3.6.3 Grain filling
A common approach is to assume grain growth as a function of biomass accumulation after anthesis plus reserves (e.g. (Jamieson et al., 1998) or simply as the sum of accumulated biomass after a defined growth stage (e.g. Ritchie 1985). An alternative approach is to simulate potential grain yield as a linear increase in harvest index (HI) with accumulated thermal time from the beginning of grain-filling to maturity (e.g. (Chapman et al., 1993; Goyne et al., 1996; Hammer et al., 1995; Meinke et al., 1998b; Sinclair and Amir, 1992) assuming a constant or cultivar specific potential HI. Thereby, the harvest index is the ratio of grain to shoot dry mass at harvest. Simplified (according to the ‘SWMIN’ concept of Ceres-Wheat) the straw fraction of the shoot at harvest can be seen as proportional to the shoot mass at begin of spike growth (W~0), where grain mass is proportional to the shoot accumulation after ear emergence ($W$):
HI=\frac{{{W}_{grain}}}{{{W}_{straw}}+{{W}_{grain}}}\sim\frac{\Delta W}{{{W}_{0}}+\Delta W} \tag{63}
The shoot mass accumulation is the product of the average amount of intercepted light (avL) times the average light use efficiency (avLUE) over the period of spike growth(\Delta t)
\Delta W=\Delta t\times avL\times avLUE \tag{64}
whereby \Delta t is the duration of the spike growth phase in thermal time (dtt) divided by the average temperature during this period (avT):
\Delta t=\frac{dtt}{avT} \tag{65}
Insertion of Eq. 63 into Eq. 65 yields:
\Delta W=\frac{dtt\cdot avL\cdot avLUE}{avT}=avQ\cdot avLUE\cdot dtt \tag{66}
whereby avQ is photo-thermal-quotient [MJ.°C-1], defined as the ratio of mean global radiation and mean temperature over the period of \Delta t.
3.6.3.1 Potential harvest index
Assuming dtt and avLUE as constant, variation in HI is a function of W0 and avQ only. In unstressed crops the variation in W0 is small compared to the variation of \Delta W, which is more influenced by water limitation. Thus, the potential HI (HIpot) is dominated by avQ. Following this reasoning, the potential HI in HumeWheat is calculated as a function of photo-thermal-quotient avQHI [MJ.°C-1] defined as the running ratio of mean global radiation and mean temperature over the period between booting and mid of milk ripe stage (BBCH 50-75):
HI_{pot}(t)=H{{I}_{\operatorname{int}}}+H{{I}_{inc}}\times av{{Q}_{HI}}(t) \tag{67}
where HI_{int} and HI_{inc} are parameters of the regression equation. The default values of HI_{int} and HI_{inc} are set to 0.196 and 0.2786 in the model.
The value of QHI is continuously calculated/updated from BBCH 60 to BBCH 75 from the ratio of sum of accumulated radiation to the sum of accumulated temperature since >= BBCH 50.
3.6.3.2 Dynamics of grain filling
According to Ritchie and Godwin (1987), the grain mass GRNWTini [mg\cdot grain^{-1}] is initialized with 3.5 at BBCH stage 65.
3.6.3.2.1 EC stages 40 to 65
In order to consider pre/early flowering growth of the ear and the spike, die grain weight is linearly increased from BBCH 40 to 65 up to the value of 3.5 mg:
\frac{GRNWT_{pl}}{dt} = \frac{GRNWT_{ini}}{1000} \cdot GPPVAR \cdot min \left( 1, \frac{EC - 40} {65 - 40} \right) - GRNWT_{pl} \tag{68}
Because the spike is assumed as part of the stem fraction, stem weight is reduced by the sum of the initial grain mass at beginning of grain filling.
\frac {dSTMWT_{pl}}{dt} = \frac {dSTMWT_{pl}}{dt} - \frac {dGRNWT_{pl}}{dt} \tag{69}
3.6.3.2.2 EC stages 65 to 75
The dynamisation of grain filling from the potential harvest index approach is done by inserting a relative harvest index or relativ potential harvest index (RHI) which is increasing during the grain filling according to a logistic growth function.
At first the parameter P5 which determines the length of the grain filling phase within the CERES model is descaled to the temperature sum for grain filling (Xm) according to:
X_m = \frac{(P5 + 21.5) } {0.05} \tag{70}
The mid of the grain filling phase is estimated from the inflection point of the logistic function of the grain filling phase using the parameter dECDP which is set to 250 [°Cd], but should be in an order of 0.5 times X~m.
X_{mid} = dECDP \tag{71}
The initial ratio of the harvest index to the potential harvest index HI/potHI (relative harvest index, RHI) is called RHIini is calculated from the ratio of the initial grain weight (GRNWTini) to the top weight of the plant (TOPWTpl) and the potential harvest index (HIpot) according to the following equation:
RHI_{ini} = \frac{ GPP \cdot GRNWT_{ini} / 1000} { TOPWT_{pl}} \cdot \frac{1}{HI_{pot}} \tag{72}
The intrinsic growth rate of the logistic increase of RHI during the grain filling phase rRHI is according to the following equation from the initial relative harvest index (RHIini) and the temperature sum of mid of the grain filling phase (Xmid)) approximated by the inflection point of the logistic growth equation:
r_{RHI} = ln \left( \frac{ 1 } {RHI_{ini} - 1 }\right) \cdot \frac{1}{Xmid} \tag{73}
The derivation of this equation is based on equation determining the time of the inflection point of a logistic growth function t_{inf} here called Xmid and the assumption that the maximum value of RHI is 1:
t_{inf}= \frac{1}{r_{RHI}} \cdot ln \left( \frac{1−RHI_{ini}}{RHI_{ini}} \right) \tag{74}
Solving equation Eq. 74 for r_{RHI} yields equation Eq. 73.
The growth rate parameter b_RHI is calculated from the potential relative harvest index (which is 1) and the initial ratio of the harvest index to the potential harvest index according to the following equation:
b_{RHI} = \frac{1 - RHI_{ini}}{ RHI_{ini}} \tag{75}
The ratio of the actual harvest index to the potential harvest index (RHI) is calculated from a logistic equation with the potential harvest index KGF and the relative growth rate of the grain filling phase according to the following equation:
RHI = \frac{1}{(1 + b_{RHI} * e^{(-r_{RHI} \cdot SUMDTT_{GF}))}} \tag{76}
See Figure 18 for the relationship between the relative potential harvest index and the cumulative temperature sum during grain filling.
plot(SUMDTTGF, RHI, type = "l", xlab = "SumDTTGF [°Cd]", ylab = "Relative harvest index (RHI)")
The grain DM growth rate GROGRN is calculated from the potential, sink limited growth rate of the grains and the difference between the actual stem weight and the minimum stem weight:
GROGRN_{pot} = (TOPWT_{pl} + Senwt_{pl}) \cdot HI_{pot} * RHI - GRNWT_{pl} \tag{77}
The growth rate of grain per plant is calculated from the potential, sink limited growth rate of the grains and the difference between the actual stem weight and the minimum stem weight SWMIN:
GROGRN = max(0, min(GROGRN_{pot}, STMWT_{pl} - SWMIN))
\frac{STMWT_{pl}}{dt} = \frac{STMWT_{pl}}{dt} - GROGRN
The grain weight per plant is calculated from the grain weight per plant and the growth rate of the grain weight per plant according to the following equation: \frac{GRNWT_{pl}}{dt} = GROGRN
4 Leaf Area Development
Many crop models are simulating leaf area as a function of simulated leaf biomass multiplied with the ratio of leaf area over leaf dry weight (SLA). Thereby two fundamentally different approaches can be employed: (i) The SLA refers to the new leaf area formed or (ii) to the overall canopy of the crop (Marcelis et al., 1998). The former approach considers LAI as a state variable, while SLA refers to the change rate of this state. With the latter approach LAI is usually defined as a derived variable. Both approaches have their justification: During early leaf growth SLA is decreasing with increasing leaf size, because of the greater proportion of structural C needed to maintain bigger leafs (Rawson et al., 1987). Thus the SLA of the new formed leafs is smaller, while the SLA of the previously formed, smaller leafs remains higher. This relationship is reflected by the decrease in SLA of new formed leaf area. As a converse effect, the influence of mutual shading becomes more important with rising LAI, what enhances SLA. The influence of shading does not concern new formed leaf area of the top leaves, but earlier formed, lower inserted leaves. Considering this aspect, it is not appropriate to conserve the SLA of already formed leaves. This implies an interaction between SLA and LAI, which should be taken into account. During early leaf growth the decrease due to leaf size probably determines SLA, but is soon dominated by the influence of mutual shading. The transition between these two phases is variable, which complicates the SLA modelling. Furthermore, for winter annual crops also the drastically changing light environment from sowing in autumn until rapid growth in spring affects SLA.
A compromising solution of both approaches is provided by the simulated leaf area development in HumeWheat. Thereby a potential SLA (SLApot) refers to the complete canopy, but affects only the LAI increase, whereas the reduction of LAI due to the senescence processes is calculated by separate routines:
LAI_{new}(t)=\max (0,W_{leaf}(t) \cdot SLA_{pot}(t) \cdot 10^{-4} - LAI(t)) \tag{78}
Daily LAI increment is calculated from the increase in total leaf area and total senescence:
\frac{dLAI}{dt}=LA{{I}_{new}}(t)-\frac{dSen}{dt} \tag{79}
4.1 Calculation of SLApot
At the beginning of leaf growth, i.e. after emergence of the first leaf SLApot is initialized with the maximum SLA (SLAmax). The standard value of SLAmax is 250 [cm2/g].
At later stages, i.e. after start of stem elongation (BBCH 31) The SLA usually increases due to mutual shading as a function of LAI.
SLA_{LAI}(t)=a_{SLA}+b_{SLA} \cdot LAI(t) \tag{80}
where a_{SLA} and b_{SLA} are parameters of the regression equation. The default values of a_{SLA} and b_{SLA} are set to 136.69 and 14.93 in the model (Figure 19). See Fig. 3b in (Ratjen and Kage, 2013).
In order to describe the transition from the high SLA values after emergence to lower values followed by the subsequent increase due to the shading effects of increasing LAI a transition value of SLA (SLATrans) is calculated:
SLA_{Trans}(t)=SLA_{LAI}(t)+ (SLA_{max}-SLA_{LAI}(t)) \cdot e^{ (f_{dec1} \cdot LAI(t)+f_{dec2})} \tag{81}
where f_{dec1} and f_{dec2} are parameters describing the transition. The default values of f_{dec1} and f_{dec2} are set to -1.1237 and 0.3 in the model.
In every time step the potential SLA is calculated as the minimum of the previous SLA and the calculated SLATrans. If the old SLA is larger than SLALAI, the SLA is equal to SLAtrans, if not SLA is equal to SLALAI.
An example for the result of this algorithm is given in Figure 20.
4.2 Leaf number
The leaf number (LN) on the main stem is calculated from the cumulative phyllochrons (CUMPH) since emergence and the phyllochron interval (PHINT) according to the following equations:
CUMPH = T_{eff} / PHINT \tag{82}
LN = trunc( CUMPH + 1 ) \tag{83}
where T_{eff} is the effective temperature, PHINT is the phyllochron interval, and LN is the leaf number on the main stem. The value of PHINT is set to 91.74 [°Cd] in the model (see Eq. 11 and Johnen et al. (2012a)).
4.3 Leaf growth
The leaf area is calculated on a leaf age class level, thereby the leaf number on the main stem is giving the age classes. This approach is basically from CERES-Wheat. Basically, only the youngest leaf age class is growing. This approach allows to implement a simple approach for leaf senescence, which is based on the leaf age. The growth rate of the leaf area of the youngest leaf class is calculated from the sum of the existing leaf dry matter of the leaf (LFWTpl) plus the daily dry matter growth rate of the leaf fraction of a single plant (GROLF) times the specific leaf area (SLA) minus the existing leaf area for that leaf age class (sumplsc):
\frac{dPLSC_i}{dt} = PLSCGR_i= (PLSC_i + GROLF) \cdot SLA - PLSC_{sum} \tag{84}
4.4 Senescence
Leaf senescence is governed by different factors in HumeWheat:
- age
- drought stress
- nitrogen stress
- light availability and also the EC stage of the crop
4.4.1 Age
4.4.1.1 CERES algorithm for leaf senescence in early stages (BBCH 10-30)
The model includes the option to simulate leaf senescence in the early stages of the crop according to the CERES-Wheat approach. This algorithm, however, was not evaluated for own data sets and the leaf area simulation may perform better if this option is set to false.
The 5th oldest leaf fraction dies within one phyllochron interval. First the leaf area of the oldest leaf fraction is fixed as a potential senescence rate when four younger leaves are present:
senratesLA[trunc(LN)-4] = PLSC[trunc(LN)-4] \tag{85}
Then the leaf loss rate of that leaf age fraction is calculated from the leaf area of this fraction at the beginning of the senescence process and the fraction of the leaf fraction of the effective temperature (Teff) and the phyllochron interval (PHINT):
PLALR_a = min(\frac { PLSC[trunc(LN) - 4]} {dt}, senratesLA[trunc(LN) - 4] * T_{eff} / (PHINT)) \tag{86}
where PLSC is the leaf area of the leaf age class, T_{eff} is the effective temperature, PHINT is the phyllochron interval, and PLALR_a is the leaf area loss rate due to age.
4.4.1.2 Leaf senescence in later stages (BBCH 31-71)
In later stages the leaf senescence is calculated as a fraction of the existing leaf area per plant, the effective day temperature and a stage specific parameter. (see https://nowlin.css.msu.edu/wheat_book/CHAPTER2.html)
For ISTAGE 2 and 4 (BBCH 30-57, see Table 1), the equation is:
PLALR_a = PSENLeaf_1 \cdot T_{eff} \cdot GPLA \tag{87}
and similarly at ISTAGE 5 (BBCH 58-71):
PLALR_a = PSENLeaf_2 \cdot T_{eff} \cdot GPLA \tag{88}
where PSENLeaf_1 is a parameter of the regression equation.
The parameters values PSENLeaf_2 and GPLA have defaults of 0.0003 and 0.0006 [cm^2 \cdot °Cd^{-1}] in the model.
4.4.2 Drought stress
The base for the calculation of the drought stress is a moving average of the relative TranspirationInterception ratio, i.e. the ratio of the sum of the actual transpiration and interception to the sum of potential transpiration and interception.
If this average ratio (TransIntRatioeven) and the actual TransIntRatio are smaller than a threshold value (TRcrit) drought stress induced leaf senescence is calculated. The threshold value is set to 0.62 in the model. This means that drought stress induced leaf senescence occurs if the average and the actual ratio of actual to potential transpiration+interception is smaller than 0.62.
In a next step a sustainable leaf area index (LAIs) is calculated which gives a ratio of the sum of actual transpiration + interception to the sum of potential transpiration to interceptionn equals the value of TRcrit.
The daily senescence rate due to drought stress then is calculated from the difference of the actual LAI and the sustainable LAI multiplied with the TransIntRatio divided by 15 (days) times the TransIntRatioeven and divided by the number of plants per m2, thereby the factor 1e4 converts m2 to cm2
PLALR_d = max\left( 0, \frac{LAI - LAIs} {15} \cdot TransIntRatio_{even} \cdot 1e4 \cdot \frac{1}{plants} \right) \tag{89}
where LAI is the leaf area index, LAIs is the sustainable leaf area index, TransIntRatio_{even} is the moving average of the relative transpiration interception ratio, and plants is the number of plants per square meter.
The sustainable leaf area index is derived from an iterative, numerical solution using the Newton method for solving non-linear equations.
4.4.3 Nitrogen stress / nitrogen induced leaf senescence
The algorithm for nitrogen stress induced leaf senescence is based on the specific leaf nitrogen concentration (SLN) and the critical specific leaf nitrogen concentration (SLNtotcrit). The daily senescence rate due to nitrogen stress is calculated from the difference of the actual LAI and the sustainable LAI multiplied with the ratio of the specific leaf nitrogen concentration to the critical specific leaf nitrogen concentration:
PLALR_n = min\left( LAImax \cdot \frac{PLALR_{max}} {100}, LAI \cdot \left( 1 - \frac{SLN} {SLNtot_{crit}} \right) \right) \cdot \frac{1e4}{plants} \tag{90}
where LAImax is the maximum leaf area index, PLALR_{max} is the maximum leaf area index loss rate [%], standard value is 5%, SLN is the specific leaf nitrogen concentration [g/m2 leaf], SLNtot_{crit} is the critical specific leaf nitrogen concentration (default 0.8) , and plants is the number of plants per square meter. The factor 1e4 converts m2 to cm2. This algorithm also induces the leaf senescence during ripening.
Alternatively, the classical leaf senecescence calculation from the CERES-Wheat model can be used. With this algorithm, the leaf area decline during ripening is controlled by the accumulated temperaturs sum during ISTAGE 5:
PLALR_n = GPLA \cdot 2 \cdot SUMDTT_5 \cdot \frac { T_{eff}}{TSum_{GF}^2} \tag{91}
where GPLA is the leaf area per plant, SUMDTT5 is the accumulated temperature sum during ISTAGE 5, T_{eff} is the temperature sum increment, and TSum_{GF}^2 is the value of unscaled parameter P5, giving the length of the grain filling period in degree days. The default value of P5 is set to 11.67 in the model. According to the equation Eq. 21 the equivalent value of TSum_{GF} is 663.4 [°Cd].
4.4.4 Light limitation
The light dependent leaf senescence is calculated from the leaf area index if the 10 day average of radiation (PARav) [MJ.m-2.d-1] is falling below a threshold value named Icrit. The threshold value is set to 1 [MJ.m-1.d-1] in the model.
Light dependent senescence only occurs between BBCH 31 and 39, the latter value is, however, a parameter describing the end of leaf growth (EClgend).
The daily senescence rate due to light limitation is calculated from the difference of the actual LAI and the sustainable LAI multiplied with the radiation sum divided by 10 days and divided by the number of plants per m2, thereby the factor 1e4 converts m2 to cm2
In a first step a sustainable leaf area index (LAIs) is calculated, which depends on the radiation sum and the critical radiation sum (Icrit) according to the following equation:
LAI_s = \frac{ln(I_{crit}) - ln(PAR_{av})} {-k_{PAR}} \tag{92}
where Icrit is the critical radiation sum, PAR_{av} is the 10 day average of photosynthetically active radiation, and k_{PAR} is the extinction coefficient for PAR radiation. The default value of kPAR is set to 0.7 in the model.
4.5 Combined senescence
In order to combine the different senescence processes the maximum of either the maximum of age ord drought senescence or the maximum of nitrogen or light dependent senescence rate is calculated. In total the senescence rate is not larger than the leaf area per plant:
PLALR = min \left( \frac{LAI*1e4}{plants}, max\left(max\left(PLALR_a,PLALR_d\right), max\left(PLALR_n, PLALR_l\right) \right)\right) \tag{93}
The algorithms for leaf area senescence are mainly taken from other models like Ceres-Wheat or I_Wheat . A stringent validation has not yet been performed.
5 Evapotranspiration
6 Soil water
7 Implementation
The model was implemented using an object oriented component library termed HUME (Kage and Stützel, 1999) on the basis of Embarcadero® Delphi® 2010 (Embarcadero Technologies, Inc., USA). The model is comprised of submodels describing plant growth and development, soil water dynamics, and evapotranspiration.
8 Miscellaneous
Temperature depended processes are calculated generally using an effective temperature (Teff) [°C], which is the mean daily air temperature [°C] (Tair) minus a base temperature (Tbase = 2 °C) (Eq. 94).
T_{eff} = \begin{cases} T_{air} - T_{base} & \text{if}~T>0\\ 0 & \text{if}~T<2 \end{cases} \tag{94}