Identification of species traits enhancing yield in wheat-faba bean intercropping: development and sensitivity analysis of a minimalist mixture model

Abstract

Aim

Cereal-legume intercropping can result in yield gains compared to monocrops. We aim to identify the combination of crop traits and management practices that confer a yield advantage in strip intercropping.

Methods

We developed a novel, parameter-sparse process-based crop growth model (Minimalist Mixture Model, M³) that can simulate strip intercrops under well-watered but nitrogen limited growth conditions. It was calibrated and validated for spring wheat (Triticum aestivum) and spring faba bean (Vicia faba) grown as monocrops and intercrops, and used to identify the most suitable trait combinations in these intercrops via sensitivity analyses.

Results

The land equivalent ratio of intercrops was greater than one over a wide range of nitrogen fertilizer levels, but transgressive overyielding, with total yield in the intercrop greater than that of either sole crop, was only obtained at intermediate nitrogen applications. We ranked the local sensitivities of the individual yields of wheat and faba bean of the whole intercrop under various nitrogen input levels to various crop traits.

Conclusions

The total intercrop yield can be improved by selecting specific traits related to phenology of both species, as well as light use efficiency of faba bean and, under high nitrogen applications, of wheat. Changes in height-related crop traits affected individual yields of species in intercrops but not the total intercrop yield.

Introduction

Intercropping is the practice of growing more than one species in the same place with a substantial overlap between their growing seasons (Willey and Rao 1980a; Willey and Rao 1980b). This is in contrast to the more common practice of growing monocrops, in which only one species is grown in a field during the growing season. Intercropping can have various advantages over growing monocrops. For instance, intercropping often results in an increased abundance of natural enemies of herbivores relative to monocrops (Letourneau et al. 2011), a stronger weed suppression in comparison with monocrops (Liebman and Dyck 1993), higher soil organic carbon and nitrogen (N) contents (Cong et al. 2015), and a higher yield stability (Raseduzzaman and Jensen 2017). Yield advantages of intercrops can occur if the growth of at least one of the species in the intercrop is less inhibited due to the combined effect of interspecific and intraspecific competition than it would be due to intraspecific competition when grown as monocrop (Klimek-Kopyra et al. 2013). This can occur if there is sufficient niche differentiation between the species in the intercrop (Jensen 1996; Li et al. 2020; Loreau 2010; Malezieux et al. 2009; Tilman 2020; Yahuza 2011).

Yield advantages of intercropping over monocrops can be substantial. One commonly used index to quantify the yield advantage of an intercrop over a monocrop is the land equivalent ratio LER (De Wit 1960; Willey and Rao 1980b). For an intercrop of two species, LER is defined as the sum of the partial land equivalent ratios of each species (pLER₁ and pLER₂, where the subscripts 1 and 2 refer to the species in the intercrop). The partial land equivalent ratio of a species is defined as the ratio of the yield of this species in an intercrop with the other species to the yield that is obtained in the monocrop, such that the LER is given by:

$$ \mathrm{LER}=p{\mathrm{LER}}_1+p{\mathrm{LER}}_2=\frac{Y_{\mathrm{ic},1}}{{\mathrm{Y}}_{\mathrm{mc},1}}+\frac{Y_{\mathrm{ic},2}}{Y_{\mathrm{mc},2}} $$

(1)

where Y_{ic, 1} and Y_{ic, 2} are the yields (weight per unit of area) of species 1 and 2 grown in an intercrop; Y_{mc, 1} and Y_{mc, 2} are the yields of species 1 and 2 obtained when they are grown as monocrops. If the LER is larger than 1, the land area that is required to obtain a certain yield with two monocrops is larger than in an intercrop. Meta-analyses report average LER values for 1,32 for maize-soybean (Xu et al. 2020), 1.29 for intercrops with maize (Li et al. 2020), 1.17 for cereal-legume intercrops (Yu et al. 2016) and, for a broader range of species, 1.30 (Martin-Guay et al. 2018). The LER is often used to quantify whether or not two species yield more in an intercrop than if these species would be grown as two separate pure cultures, but a LER > 1 does not necessarily indicate that the total yield of the intercrop is higher than the highest yielding monocrop; i.e., that transgressive overyielding (Cardinale et al. 2007; Fridley 2001; Loreau 2010). To quantify whether the intercrop confers an advantage in terms of total yield relative to a monocrop of the highest yielding species in that that intercrop, the transgressive overyielding index (TOI) should be used (Cardinale et al. 2007; Yu 2016):

$$ \mathrm{TOI}=\frac{Y_{\mathrm{ic},1}+{Y}_{\mathrm{ic},2}}{\max \left({Y}_{\mathrm{mc},1},{Y}_{\mathrm{mc},2}\right)} $$

(2)

where the max() function takes the maximum of yields of both species that would be obtained in a monocrop, and, as above, Y is the yield, in intercrop (subscript ic) and monocrop (subscript mc), for the two species (1 and 2).

Cereal-legume intercrops are an example of an intercropping system with a clear potential for niche differentiation. While cereals take up inorganic nitrogen (NO₃⁻, NH₄⁺) from the soil to fulfil their nitrogen demands, legumes fulfill a large part of their nitrogen demand by fixing atmospheric nitrogen (N₂). Consequently, the cereal plants experience less competition for inorganic nitrogen from a legume neighbor than from a neighbor of its own species. Yet, how advantageous this reduced competition for nitrogen is, depends on nitrogen availability in the soil, and hence the nitrogen input level. Furthermore, cereal and legume species in an intercrop still compete for other resources like water (Ozier-Lafontaine et al. 1998), solar radiation (Keating and Carberry 1993) and nutrients other than nitrogen. This remaining competition could substantially reduce the advantages of complementarity for nitrogen acquisition of growing these species together. In an intercrop, the extent of interspecific competition for these other resources is determined by differences in various crop traits between the two species. For example, differences in canopy structure and crop height (Gou et al. 2017a; Gou et al. 2017b; Keating and Carberry 1993; Kropff 1993; Pronk et al. 2003), root system architecture (Corre-Hellou et al. 2007; Ozier-Lafontaine et al. 1998), and nutrient uptake capacity (Corre-Hellou et al. 2006) substantially affect the competitive ability of one species against the other one. Additionally, management decisions like row configuration, plant density, fertilization amounts and timing, as well as sowing dates and densities can also play a role in the interspecific competition (Yu et al. 2015). All these variables affect interspecific competition simultaneously and, therefore, yield, TOI, and LER. Hence, a system approach is essential to study a cereal-legume system.

Crop growth models allow a system approach to the evaluation of intercropping systems. These process-based models can be parametrized and applied to a variety of conditions and thus allow evaluating the performance of cereal and legume species with a wide range of traits in monocrops and intercrops grown under different pedo-climatic conditions and nitrogen availabilities. The extensive analysis of the combined effects of crop traits and pedo-climatic conditions made possible by a model are a stepping-stone towards the development of varieties better suited for applications in intercrops. Nevertheless, only few crop models have been developed to simulate intercropping systems. Examples are APSIM (Keating et al. 2003), CERES (Hoogenboom et al. 2017; Ritchie and Otter 1985) and STICS (Brisson et al. 2003). Even for these models, applications to intercrops are few. For some examples, see Brisson et al. (2004) and Corre-Hellou et al. (2009) for STICS; Carberry et al. (1996), Chimonyo et al. (2016), Knörzer et al. (2011b), Weih et al. (2019), and Berghuijs et al. (2019) for APSIM; Knörzer et al. (2011b) and Knörzer et al. (2011a) for the model CERES. In the case of APSIM, the model was found to have limited capabilities to simulate cereal-legume intercrops of wheat-faba bean (Berghuijs et al. 2019; Weih et al. 2019) and wheat-fieldpea (Knörzer et al. 2011b), because the competitive ability of one species in an intercrop over the other was overestimated. Furthermore, the aforementioned models exclusively consider homogeneously mixed system and do not account for the heterogeneity of light capture in strip intercropping systems. A strip intercrop is defined as an intercrop, in which the species are grown in separate strips consisting of one or more crop rows.

Aside for the geometric arrangement of the plants, simulating intercrops with these models can be challenging, because they contain a large number of parameters. The uncertainties in the estimates of these parameters are large given the set of measurements typically available from most field experiments. This makes it difficult to use them to simulate existing intercropping systems and determine which combinations of, for instance, a cereal and a legume cultivar result in yield advantages.

A minimalist modelling approach can reduce the uncertainties inherent in parameter estimation. Moreover, a limited number of parameters makes it easier to adjust the model for other combinations of species in an intercrop grown under different conditions (Van der Werf et al. 2007). Some minimalist models have been developed for strip intercrops of wheat-maize under potential growth conditions (Gou et al. 2017b; Liu et al. 2017) and water limited growth conditions (Tan et al. 2020). Nevertheless, to the best of our knowledge, there are currently no crop growth models for any species to simulate strip intercrops grown under nitrogen-limited growth conditions. In this study, we developed, calibrated, and validated a process-based minimalist model (M³: Minimal Mixture Model). We aimed to address the need for a model that is suitable for strip intercropping and accounts for nitrogen limitation, while requiring as few parameters as possible, while keeping the model suitable for making simulations that are useful for making predictions under actual field conditions. M³ explicitly describes intra-and inter-specific interactions typical of strip intercrops and relies on a limited number of parameters, which have clear bio-physical meanings. It aims to simulate cereal-legume intercrops explicitly considering, amongst others, nitrogen limited growth and row configuration. The model is used to answer the following research questions, of relevance for management and variety choice and breeding: 1) for which nitrogen input levels is the total yield of a cereal-legume strip intercrop higher than the monocrop yield of any of the species present in this intercrop? 2) which crop traits affect the total grain yield of a strip cereal-legume intercrop the most? 3) which crop traits affect the individual yields of the cereal and the legume in a cereal-legume strip intercrop the most? 4) how does the nitrogen input level influence which of the crop traits affect grain yields the most? In the remaining sections, we first describe the development of M³ and its calibration and validation based on data of monocrops and intercrops of spring wheat (Triticum aestivum) and spring faba bean (Vicia faba) grown in Wageningen (The Netherlands). Finally, we present the results of local and global sensitivity analyses to answer the previously stated research questions.

Materials and methods

Model description

General overview M³

M³ simulates the daily biomass production and its partitioning over plant organs of either one crop species (monocrop) or two crop species (intercrop) that have partially or fully overlapping growing seasons. The state variables of the model are i) the leaf area index (i.e., the projected leaf are per unit ground area; L_i(t); m² leaf m^/2 ground), ii) the total aboveground biomass (B_i(t); kg dry matter m⁻² ground), iii) the storage organ dry matter (Y_i(t); kg dry storage organ matter m⁻² ground), iv) crop height (H_i(t); m), v) the temperature sum from sowing (S_i(t); ^°C d), vi) the crop nitrogen amount (N_{c, i}(t); kg N m⁻² ground), and vii) the soil mineral nitrogen amount in the rooting zone (N_s(t); kg mineral N m⁻² ground). The subscript i refers to the crop species. In the remainder of this manuscript, when necessary, it is replaced by “w” or “f”, in order to indicate whether it refers to wheat or faba bean.

M³ considers the joint effects of two potential limitations to plant growth: photosynthetically active radiation and nitrogen, explicitly describing the light interception by the crops and the soil and plant nitrogen dynamics, and their effects on crop growth. The focus on nitrogen stems from the pivotal role of nitrogen availability for plant growth and the niche differentiation typical of a legume-cereal intercrop. Conversely, M³ does not simulate the effect of water stress on plant growth, thus implicitly assuming that the plants are well-watered throughout the season in which drought stress was prevented. In line with this assumption, all experiments that were used to calibrate and validate M³ were irrigated. Plants extract resources, nitrogen in particular, from the rooting zone i.e. the soil layer where most of the roots are located. The depth of this zone is assumed the same for both crop species, and characterized by horizontally and vertically uniform conditions. Given the lack of information on root distribution in most field experiments, this depth is assumed constant throughout the growing season.

M³ calculates daily how much photosynthetic active radiation is intercepted by either of the two species in an intercrop or by a single species in a monocrop. From that, M³ calculates the potential daily biomass production (i.e., the production in the absence of water or nutrient shortages, pest and diseases or other growth reducing factors) (Van Ittersum and Rabbinge 1997) of each species.

Subsequently, M³ calculates for each species the demand for inorganic nitrogen. This demand for nitrogen is the amount of nitrogen that needs to be taken up by the plant in order to sustain a nitrogen concentration in the tissues that would be more than enough to enable maximum performance. This concentration is referred to here as the maximum nitrogen concentration of the crop (Justes et al. 1994; Stöckle and Debaeke 1997); see Fig. 1 for a schematic overview. Legumes can obtain a fraction of their demand by nitrogen fixation. If there is enough soil mineral nitrogen in the soil to meet the demand of a species, that species takes up an amount of nitrogen equal to that demand. If its demand cannot be met, the species takes up all nitrogen it has access to. If the crop nitrogen concentration in the plant drops below a critical crop nitrogen concentration (Greenwood et al. 1990), the actual daily biomass production is reduced relative to the potential daily biomass production (Fig. 2). Note that the maximum nitrogen concentration is larger than the critical nitrogen concentration, which allows M³ to consider that crops can take up more nitrogen then they need to reach their maximum relative growth rate (Greenwood et al. 1990; Justes et al. 1994; Van Wijk et al. 2003; Ågren 1988; Ågren and Weih 2020). If the crop nitrogen concentration drops below a so-called minimum nitrogen concentration, the actual daily biomass production is reduced to zero. Even if the plants are grown under conditions in which they can fulfill their nitrogen demand throughout their growing season, the crop nitrogen concentration still decreases over time, because the leaf:stem weight ratio decreases over time and stems maintain lower nitrogen concentrations than leaves (Greenwood et al. 1990; Justes et al. 1994). M³ considers this nitrogen dilution by assuming empirical power relationships between the aboveground dry matter and the minimum, critical and maximum nitrogen concentrations; these relationships are often referred to as “dilution curves” (Justes et al. 1994). After calculation of the newly produced biomass for a species, the newly created biomass is distributed over the leaves, storage organs, and other plant organs of that species.

Fig. 1

Schematic representation of nitrogen dilution curves (Greenwood et al. 1990; Justes et al. 1994; Stöckle and Debaeke 1997) of the maximum (ν_{max, i}(t)), critical (ν_{crit, i}(t)) and minimum nitrogen concentrations (ν_{min, i}(t)) with aboveground dry crop biomass. At very low biomasses, the maximum, critical and minimum nitrogen concentrations are equal to ν_{max, 0, i}, ν_{max, 0, i} ∙ f_{crit, i}, and ν_{max, 0, i} ∙ f_{min, i} respectively. At higher biomasses, the curves are described by a power law. The nitrogen demand D_i(t) for a given biomass B_i(t) is calculated as the difference between the amount of nitrogen that would result to the maximum nitrogen concentration for that biomass, ν_{max, i}(t) ∙ B_i(t)

Fig. 2

Schematic representation of the relationship between the growth reduction factor due to nitrogen stress (σ_{N, i}(t)) and the crop nitrogen concentration. The growth reduction factor is 0 if the crop nitrogen concentration is smaller than the minimum nitrogen concentration ν_{min, i}(t). For nitrogen concentrations between the minimum and the critical nitrogen concentration ν_{crit, i}(t), the growth reduction factor increases from 0 to 1. If the crop nitrogen concentration is higher than ν_{crit, i}, the nitrogen growth reduction factor remains constant at 1

Finally, M³ calculates daily biomass loss, and the accompanying nitrogen loss, due to leaf senescence. Re-allocation is considered by assuming the nitrogen concentration of dying leaves equals the minimum nitrogen concentration. In parallel, M³ defines the temporal evolution of phenological stage, leaf area index and plant height.

The equations describing the temporal evolution of the state variables are justified in the next sections. The model was implemented in the C# programming language (Hejlsberg et al. 2010) within the .NET 4.6.1 framework (Microsoft Corporate). The dynamics of the key variables is forced with observed meteorological data and solved numerically using an explicit Euler scheme with a daily time step (i.e., Δt = 1 d). The source code will be published on the website of the EU Horizon 2020 project Diversify (https://www.plant-teams.eu/).

The remaining part of this section describes M³ in detail.

Simulation of developmental stage

Temperature accumulation drives phenological phenomena like dry matter partitioning, flowering, and grain filling. Hence, temperature sum is calculated. The rate of increase of the temperature sum for crop i, S_i, from the date of sowing is calculated as (Van der Werf et al. 2007):

$$ \frac{\mathrm{d}{S}_{\mathrm{i}}(t)}{\mathrm{d}t}=\max \left(0,T(t)-{T}_{\mathrm{b},\mathrm{i}}\right) $$

(3)

where T(t) is the temperature (°C) at time t; and T_{b, i} is the base temperature (°C) of species i, i.e., the threshold temperature below which no increase in temperature sum occurs. The temperature sum from emergence is used to calculate the development stage δ_i(t) as:

$$ {\delta}_{\mathrm{i}}(t)=\left\{\begin{array}{ccc}0& \mid & {S}_{\mathrm{i}}(t)\le {S}_{\mathrm{e},\mathrm{i}}\\ {}\frac{S_{\mathrm{i}}(t)-{S}_{\mathrm{e},\mathrm{i}}}{S_{\mathrm{a},\mathrm{i}}}& \mid & {S}_{\mathrm{e},\mathrm{i}}<{S}_{\mathrm{i}}(t)\le {S}_{\mathrm{e},\mathrm{i}}+{S}_{\mathrm{a},\mathrm{i}}\\ {}1+\frac{S_{\mathrm{i}}(t)-{S}_{\mathrm{a},\mathrm{i}}-{S}_{\mathrm{e},\mathrm{i}}}{S_{\mathrm{m},\mathrm{i}}}& \mid & {S}_{\mathrm{e},\mathrm{i}}+{S}_{\mathrm{a},\mathrm{i}}<{S}_{\mathrm{i}}(t)\le {S}_{\mathrm{e},\mathrm{i}}+{S}_{\mathrm{a},\mathrm{i}}+{S}_{\mathrm{m},\mathrm{i}}\\ {}2& \mid & {S}_{\mathrm{i}}(t)\ge {S}_{\mathrm{e},\mathrm{i}}+{S}_{\mathrm{m},\mathrm{i}}+{S}_{\mathrm{a},\mathrm{i}}\end{array}\right. $$

(4)

where S_{e, i}, S_{a, i} and S_{m, i} are the temperature sums required to transition from sowing to emergence, emergence to anthesis and anthesis to maturity, respectively. δ_i(t) ranges between 0 and 2: it is 0 until emergence, 1 at anthesis and 2 after physiological maturity (Supit and Hooijer 1994). Anthesis refers to 50% flowering in both wheat and faba bean. Figure 3 shows eq. 4 graphically.

Fig. 3

Schematic overview of the relationship between development stage and temperature sum from sowing (^°C d)

Simulation of leaf area index

The leaf area index of species i, L_i(t), changes in time as the net result of leaf biomass growth, driven by light and nitrogen availability, and senescence:

$$ \frac{\mathrm{d}{L}_{\mathrm{i}}(t)}{\mathrm{d}t}={\sigma}_{\mathrm{N},\mathrm{i}}(t)\cdotp {S}_{\mathrm{la},\mathrm{i}}\cdotp {f}_{\mathrm{L},\mathrm{i}}(t)\cdotp {P}_{\mathrm{pot},\mathrm{i}}(t)-{\mu}_{\mathrm{i}}(t)\cdotp {L}_{\mathrm{i}}(t) $$

(5)

where σ_{N, i}(t) is a factor between 0 and 1 that accounts for a reduction of crop growth of species i at time t as a result of nitrogen limitation (i.e., if the crop nitrogen concentration is below its critical nitrogen concentration), which will be further specified. S_{la, i} is the specific leaf area (m² leaf kg⁻¹ leaf) of species i. f_{L, i}(t) is the fraction of newly produced biomass of species i allocated to the leaves at time t. P_{pot, i}(t) is the rate of total dry matter production of species i (kg m⁻² d⁻¹), which is calculated as:

$$ {P}_{\mathrm{p}\mathrm{ot},\mathrm{i}}(t)={\eta}_{\mathrm{i}}\cdotp {f}_{\mathrm{p}}\cdotp {f}_{\operatorname{int},\mathrm{i}}(t)\cdotp {I}_0(t) $$

(6)

where η_i is the light use efficiency of photosynthetic active radiation (PAR) by species i. f_p is the fraction of photosynthetic active radiation (PAR) in the global radiation (J J⁻¹, assumed to be 0.5 (Van Oijen and Leffelaar 2008)). f_{int, i}(t) is the fraction of global radiation that is intercepted by species i at time t. I₀(t) is the global radiation (J radiation m⁻² land d⁻¹) at time t. μ_i(t) is the relative mortality rate of leaf biomass (d⁻¹). The calculation of f_{int, i}(t) largely follows the radiation interception model described by Gou et al. (2017a). This model takes into account the light distribution in a mixed canopy with species grown in strips, whereby the strip width, the plant height and the leaf area index per species are accounted for using algorithms first derived by Goudriaan (1977) and Pronk et al. (2003). A detailed description can be found in Supplementary Text S1.

The calculation of the biomass partitioning to leaves in M³, f_{L, i}(t), follows the approach of the Yield-SAFE model (Van der Werf et al. 2007). Until species i reaches the developmental stage δ_i(t) = d_{1, i} (see equation 4), f_{L, i}(t) has a constant value f_{L, 0, i}. After this developmental stage has been reached, f_{L, i}(t) decreases linearly with development until it reaches the value 0 at developmental stage δ_i(t) = d_{2, i}. Beyond this development stage, no more newly produced biomass is allocated to the leaves. f_{L, i}(t) is calculated as:

$$ {f}_{\mathrm{L},\mathrm{i}}(t)=\left\{\begin{array}{ccc}{f}_{\mathrm{L},0,\mathrm{i}}& \mid & {\delta}_{\mathrm{i}}(t)\le {d}_{1,\mathrm{i}}\\ {}{f}_{\mathrm{L},0,\mathrm{i}}\cdot \frac{d_{2,\mathrm{i}}-{\delta}_{\mathrm{i}}(t)}{d_{2,\mathrm{i}}-{d}_{1,\mathrm{i}}}& \mid & {d}_{1,\mathrm{i}}<{\delta}_{\mathrm{i}}(t)\le {d}_{2,\mathrm{i}}\\ {}0& \mid & {\delta}_{\mathrm{i}}(t)>{d}_{2,\mathrm{i}}\end{array},\operatorname{}\right. $$

(7)

For the relative mortality rate μ_i(t), M³ assumes that senescence in species i takes place after anthesis (δ_i(t) > 1) and above a threshold temperature T_{sen, i}. Above this temperature, the rate of mortality due to leaf senescence increases with temperature:

$$ {\mu}_{\mathrm{i}}(t)=\left\{\begin{array}{ccc}0& \mid & {\delta}_{\mathrm{i}}(t)\le 1\vee T(t)\le {T}_{\mathrm{sen},\mathrm{i}}\\ {}\min \left(1,{s}_{\mathrm{sen},\mathrm{i}}\cdotp \left(T(t)-{T}_{\mathrm{sen},\mathrm{i}}\right)\right)& \mid & {\delta}_{\mathrm{i}}(t)>1\wedge T(t)>{T}_{\mathrm{sen},\mathrm{i}}\end{array}\right. $$

(8)

Simulation of the aboveground biomass

The rate of change of the aboveground dry matter is calculated as:

$$ \frac{\mathrm{d}{B}_{\mathrm{i}}(t)}{\mathrm{d}t}={\sigma}_{\mathrm{N},\mathrm{i}}(t)\cdotp {P}_{\mathrm{pot},\mathrm{i}}(t)-{\mu}_{\mathrm{i}}(t)\cdotp \frac{L_{\mathrm{i}}(t)}{S_{\mathrm{la},\mathrm{i}}} $$

(9)

which combines the production from equation 5 and the loss from equation 8. M³ assumes that only leaf tissue can die within the growing season: hence, the mortality rate of leaf dry matter equals the mortality rate of total aboveground dry matter.

Simulation of storage organ weights

The growth rate of the dry matter of the storage organs of species i (Y_i(t)) is calculated as:

$$ \frac{\mathrm{d}{Y}_{\mathrm{i}}(t)}{\mathrm{d}t}={\sigma}_{\mathrm{N},\mathrm{i}}(t)\cdotp {f}_{\mathrm{y},\mathrm{i}}(t)\cdotp {P}_{\mathrm{pot}}(t) $$

(10)

where f_{y, i}(t) is the fraction of newly produced biomass that is partitioned to the storage organs of species i at time t. Until a development stage δ_i(t) = d_{3, i} (equation 4) has been reached, no biomass is partitioned to the storage organs. d_{3, w} refers to the development state at which floral initiation takes place. d_{3, f} refers to the development state where the first flowers are initiated. Subsequently, the fraction of newly produced biomass that is assigned to the storage organs increases linearly with δ_i(t) until the development stage δ_i(t) = d_{4, i} is reached. After that, all newly produced biomass is partitioned to the grains and, if d_{4, i} < d_{2, i} (and consequently f_{L, i}(t) > 0), also to the leaves (Fig. 4). In summary,

Fig. 4

Schematic representation of the relationship between changes in partitioning fractions to the leaves (“Leaf”), storage organs (“Stor.”), and remaining aboveground organs (“Rem.”) and development state (equation 7 and 11). f_{L, 0} is the partitioning fraction of newly produced biomass to the leaves from emergence until δ_i(t) = δ_{1, i} As there is no biomass assigned to the storage organs at that stage, the fraction of biomass assigned to the remaining organs at emergence is 1 − f_{L, 0}

$$ {f}_{\mathrm{y},\mathrm{i}}(t)=\left\{\begin{array}{ccc}0& \mid & {\delta}_{\mathrm{i}}(t)<{d}_{3,\mathrm{i}}\\ {}\left(1-{f}_{\mathrm{L},\mathrm{i}}(t)\right)\cdotp \frac{\delta_{\mathrm{i}}(t)-{d}_{3,\mathrm{i}}}{d_{4,\mathrm{i}}-{d}_{3,\mathrm{i}}}& \mid & {d}_{3,i}<{\delta}_{\mathrm{i}}(t)<{d}_{4,\mathrm{i}}\\ {}1-{f}_{\mathrm{L},\mathrm{i}}(t)& \mid & {\delta}_{\mathrm{i}}(t)>{d}_{4,\mathrm{i}}\end{array}\right. $$

(11)

Simulation of crop height

M³ follows the model of Gou et al. (2017b) by assuming that the height growth rate of species i is proportional to the effective temperature (difference between T(t) and T_{b, i}) and that the height asymptotically approaches a maximum height H_{m, i} following a logistic growth trajectory in thermal time:

$$ \frac{\mathrm{d}{H}_{\mathrm{i}}(t)}{\mathrm{d}t}={r}_{\mathrm{h},\mathrm{i}}\cdotp \left(T(t)-{T}_{\mathrm{b},\mathrm{i}}\right)\cdotp {H}_{\mathrm{i}}(t)\cdotp \left(1-\frac{H_{\mathrm{i}}(t)}{H_{\mathrm{m},\mathrm{i}}}\right) $$

(12)

where r_{h, i} is the temperature sensitive relative growth rate of height growth (^°C⁻¹ d⁻¹) of species i.

Simulation of soil and crop mineral nitrogen

Nitrogen dynamics of the soil-plant system is a complex process. In order to keep the number of parameters low, M³ adopts a number of simplifying assumptions similar to the LINTUL3 model (Shibu et al. 2010). M³ assumes that all inorganic nitrogen that is applied as fertilizer is available for uptake by each crop species and that the net mineralization rate is constant throughout the growing season. In reality, mineralization rate is temperature dependent (Kirschbaum 1995; Kätterer et al. 1998). However, no large seasonal change in mineralization rate are expected, because we only simulate spring crops. So, the range of possible daily temperatures is limited (Eckersten et al. 2010). M³ then calculates the rate of change of the soil mineral nitrogen as the difference between the two source terms, the addition of nitrogen to the soil by fertilization application and the net mineralization, and one loss term, the removal of soil mineral nitrogen by crop uptake:

$$ \frac{\mathrm{d}{N}_{\mathrm{s}}(t)}{\mathrm{d}t}={f}_{\mathrm{recov}}\cdot F(t)+M-\sum \limits_j{U}_{\mathrm{j}}(t) $$

(13)

where f_recov is the recovery fraction, i.e., the fraction of mineral nitrogen in the fertilizer that is not lost due to volatilization and leaching; it is set to 0.7 (Wolf 2012). F(t) is the amount of mineral nitrogen that is applied at time t (kg N m⁻² d⁻¹). M is the mineralization rate (kg N m⁻² d⁻¹); we determined its value based on experimental data as the ratio of an assumed amount 25 kg N ha⁻¹ that was mineralized during the growing season (Gou et al. 2017c) and the mean of the length of the spring wheat growing season (148 d), leading to M = 1.69 · 10⁻⁵ kg N m⁻² d⁻¹. $\sum \limits_j{U}_{\mathrm{j}}(t)$ is the sum of the nitrogen uptake rates of all species in the intercrop system (kg N m⁻² d⁻¹). The calculation of uptake rates will be specified in the next section.

The change in crop nitrogen amount is simulated as:

$$ \frac{\mathrm{d}{N}_{\mathrm{c},\mathrm{i}}(t)}{\mathrm{d}t}={U}_{\mathrm{i}}(t)+{\Phi}_{\mathrm{i}}(t)-{\mu}_{\mathrm{N},\mathrm{i}}(t) $$

(14)

where U_i(t) is the uptake rate of species i (kg N m⁻² d⁻¹). Φ_i(t) is the amount of N₂ fixation at time t (kg N m⁻² d⁻¹). Both U_i(t) and Φ_i(t) will be further specified. μ_N(t) is the nitrogen loss rate (kg N m⁻² d⁻¹) of species i at time t, which will also be further specified.

Calculation of the plant nitrogen uptake and use

Four key nitrogen crop concentrations are used in M³ (Fig. 1): they are all expressed in kg N kg⁻¹ aboveground dry weight. The actual nitrogen concentration ν_i(t) is the crop nitrogen concentration at time t. The maximum nitrogen concentration ν_{max, i}(t) is the crop nitrogen concentration above which species i does not take up any nitrogen from the soil at time t. The critical nitrogen concentration is the crop nitrogen concentration ν_{crit, i}(t) above which the plant can maintain its potential growth (i.e., σ_N = 1; see equations 5, 9, and 10). If ν_i(t) is lower than ν_{crit, i}(t), the actual biomass production is smaller than the potential biomass production (σ_N(t) < 1). If ν_i(t) is also lower than the minimum nitrogen concentration ν_{min, i}(t), the biomass growth is completely halted (σ_N = 0) (Fig. 2). M³ calculates the crop nitrogen concentration ν_i(t) as:

$$ {\nu}_{\mathrm{i}}(t)=\frac{N_{\mathrm{c},\mathrm{i}}(t)}{B_{\mathrm{i}}(t)} $$

(15)

In order to calculate ν_{max, i}(t), ν_{crit, i}(t), and ν_{min, i}(t), M³ applies three nitrogen dilution curves. Previous dilution curves (Greenwood et al. 1990; Justes et al. 1994; Stöckle and Debaeke 1997) were defined for monocrops. However, in strip intercrops each species is grown on a fraction of the land. Consequently, the aboveground dry matter term in nitrogen dilution curves for a species in a strip intercrop should be expressed in aboveground dry weight per unit of land on which this species is grown. This variable is obtained from the biomasses in the model, which are expressed in kg per unit area of the system as a whole, by multiplying the biomass of that species with the inverse of the fraction of land on which this species is grown. This results in the following nitrogen dilution curves:

$$ {\nu}_{\max, \mathrm{i}}(t)={\nu}_{\max, 0,\mathrm{i}}\cdotp \min \left(1,{\left(\frac{\ell_{\mathrm{i}}+{\ell}_{\mathrm{j}}}{\ell_{\mathrm{i}}}\cdotp \frac{B_{\mathrm{i}}(t)}{B_{\max, \mathrm{i}}}\right)}^{-{\beta}_{\max, \mathrm{i}}}\ \right) $$

(16)

$$ {\nu}_{\mathrm{crit},\mathrm{i}}(t)={f}_{\mathrm{crit},\mathrm{i}}\cdotp {\nu}_{\max, 0,\mathrm{i}}\cdotp \min \left(1,{\left(\frac{\ell_{\mathrm{i}}+{\ell}_{\mathrm{j}}}{\ell_{\mathrm{i}}}\cdotp \frac{B_{\mathrm{i}}(t)}{B_{\mathrm{crit},\mathrm{i}}}\right)}^{-{\beta}_{\mathrm{crit},\mathrm{i}}}\ \right) $$

(17)

$$ {\nu}_{\min, \mathrm{i}}(t)={f}_{\min, \mathrm{i}}\cdotp {\nu}_{\max, 0,\mathrm{i}}\cdotp \min \left(1,{\left(\frac{\ell_{\mathrm{i}}+{\ell}_{\mathrm{j}}}{\ell_{\mathrm{i}}}\cdotp \frac{B_{\mathrm{i}}(t)}{B_{\min, \mathrm{i}}}\right)}^{-{\beta}_{\min, \mathrm{i}}}\ \right) $$

(18)

where are the β_{crit, i}, β_{max, i}, and β_{min, i} are dimensionless shape exponents. ν_{max, 0, i} (kg N kg⁻¹ aboveground dry matter) is the maximum crop nitrogen concentration at low biomass (in this study, we assumed a biomass less than 1000 kg ha⁻¹) for which the maximum nitrogen concentration is biomass independent (Justes et al. 1994; Stöckle and Debaeke 1997). f_{crit, i} and f_{min, i} (both <1) are the ratios of the critical and minimum concentration at these biomasses to ν_{max, 0} respectively. The parameters B_{max, i}, B_{crit, i}, and B_{min, i} are the aboveground biomasses for which, respectively, the maximum, critical and minimum nitrogen contents equal to 1, in case the model would not assume constants values for these nitrogen contents at low biomasses (i.e. ν_{max, 0, i}, f_crit · ν_{max, 0, i}, and f_min · ν_{max, 0, i}). ℓ_i is the strip width of species i and ℓ_j is the strip width of the second species in the intercrop. In monocrop of species i, ℓ_j = 0 and, therefore, the ratio $\frac{\ell_{\mathrm{i}}+{\ell}_{\mathrm{j}}}{\ell_{\mathrm{i}}}$ equals 1. Figure 1 shows the nitrogen dilution curves schematically.

The nitrogen demand D_i(t) (kg N m⁻² d⁻¹) of each species i in the system is calculated as:

$$ {D}_{\mathrm{i}}(t)={t}_{\mathrm{c}}\cdotp \max \left(0,{B}_i(t)\cdotp \left({\nu}_{\max, \mathrm{i}}(t)-{\nu}_i(t)\right)\right) $$

(19)

where t_c is a relative rate (d⁻¹) equal to 1 d⁻¹.

Whether the crop nitrogen demand can be met depends on the soil mineral nitrogen supply, the nitrogen demand of all species in the system, and their capability to fix nitrogen. M³ assumes that each species can fix a fraction f_{fix, i} of its demand from N₂ in the atmosphere (Wolf 2012). M³ calculates rate of nitrogen fixation as:

$$ {\Phi}_{\mathrm{i}}(t)={f}_{\mathrm{fix},\mathrm{i}}\cdotp {D}_{\mathrm{i}}(t) $$

(20)

f_{fix, i} equals 0 for species unable to fix N along this pathway. The extent to which faba bean obtains nitrogen from fixation relative to obtaining nitrogen from the soil varies substantially. A meta analysis of 18 studies conducted in different parts in the world (Jensen et al. 2010) showed that the fraction of shoot nitrogen obtained from fixation varies from 0.34 to 0.88 in fertilized field trials and from 0.44 to 0.99 in unfertilized field trials. It has also been shown in pea that the fraction of fixed nitrogen can vary in time (Corre-Hellou et al. 2006; Liu et al. 2011) and that it can be both positively (Ghaley et al. 2005) and negatively affected by fertilizer input (Corre-Hellou et al. 2006). Given these large uncertainties, we decided to assume f_{fix, f} to be a constant and set it equal to 0.8 (Wolf 2012), which falls within the range of values reported by Jensen et al. (2010).

All species have to fulfil the remaining part of their demand nitrogen uptake from the soil. This nitrogen demand from the soil D_{soil, i}(t) is calculated as:

$$ {D}_{\mathrm{soil},\mathrm{i}}(t)=\left(1-{f}_{\mathrm{fix},\mathrm{i}}\right)\cdotp {D}_{\mathrm{i}}(t) $$

(21)

M³ daily splits the soil mineral nitrogen into two pools for species i and j, which both contain a fraction of the total amount of soil mineral nitrogen (f_{acN, i}(t) and f_{acN, j}(t)). These fractions are determined by the transpiration rate of species i relative to species j. The transpiration rate τ(t) (m³ H₂O m⁻² soil) of species i is calculated as (Allen et al. 1998; Tan et al. 2020):

$$ {\tau}_{\mathrm{i}}(t)={\sigma}_{\mathrm{w},\mathrm{i}}(t)\cdotp {K}_{\mathrm{tr},\mathrm{i}}(t)\cdotp {E}_0(t)\cdotp {f}_{\operatorname{int},\mathrm{i}}(t) $$

(22)

where σ_{w, i}(t) is the fraction of reduction of growth and transpiration due to water stress of species i. K_{tr, i}(t) is the potential transpiration coefficient of species i at time t; this is the ratio of the potential transpiration rate of crop i to a reference crop as defined by the FAO (Allen et al. 1998). E₀(t) is the reference evapotranspiration rate (m³ H₂O m⁻² ground d⁻¹) at time t. It is assumed that σ_{w, i}(t) = σ_{w, j}(t) = 1 (no water stress, in line with the experimental setup used for model calibration and validation) and that K_{tr, i}(t) = K_{tr, j}(t). Given these assumptions, f_{acN, i}(t) is given by:

$$ {f}_{\mathrm{acN},\mathrm{i}}(t)=\frac{\tau_{\mathrm{i}}(t)}{\tau_{\mathrm{i}}(t)+{\tau}_{\mathrm{j}}(t)}=\frac{f_{\operatorname{int},\mathrm{i}}(t)}{f_{\operatorname{int},\mathrm{i}}(t)+{f}_{\operatorname{int},\mathrm{j}}(t)} $$

(23)

If species i is grown in a monocrop, f_{int, j}(t) = 0; otherwise f_{int, j}(t) is the fraction of light interception by the second species in the intercrop.

The nitrogen uptake of species i equals its nitrogen demand from the soil, if that is lower than the soil mineral nitrogen that is available to this species. If this is not the case, each species takes up all soil mineral nitrogen that they have access to:

$$ {U}_{\mathrm{i}}(t)=\left\{\begin{array}{ccc}{D}_{\mathrm{s}\mathrm{oil},\mathrm{i}}(t)& \mid & {D}_{\mathrm{s}\mathrm{oil},\mathrm{i}}(t)\le {f}_{\mathrm{acN},\mathrm{i}}(t)\cdotp {N}_{\mathrm{s}}(t)\\ {}{f}_{\mathrm{acN},\mathrm{i}}(t)\cdotp {N}_{\mathrm{s}}(i)& \mid & {D}_{\mathrm{s}\mathrm{oil},\mathrm{i}}(t)>{f}_{\mathrm{acN},\mathrm{i}}(t)\cdotp {N}_{\mathrm{s}}(t)\end{array}\right. $$

(24)

Each day the reduction factor of biomass production is calculated as a piecewise linear function of the plant nitrogen concentration at that time:

$$ {\sigma}_{\mathrm{N}}(t)=\max \left(0,\min \left(1,\frac{\nu_i(t)-{\nu}_{\min, \mathrm{i}}(t)}{\nu_{\mathrm{crit},\mathrm{i}}(t)-{\nu}_{\min, \mathrm{i}}(t)}\right)\right) $$

(25)

where biomass production is completely inhibited if σ_N(t) = 0, i.e., if ν_i(t) ≤ ν_{min, i}(t). Conversely, the biomass production is equal to the potential production rate if σ_N(t) = 1 (Fig. 2).

Crop nitrogen loss

M³ assumes that crops lose nitrogen via leaf senescence while the senescence of other plant organs is negligible. M³ considers that the plants minimize this nitrogen loss by nitrogen retranslocation from senescing leaves to other plant organs, so that the nitrogen concentration of senesced leaves is ν_{min, i}(t). As such, M³ calculates the rate of nitrogen loss μ_N (see equation 14) as:

$$ {\mu}_{\mathrm{N}}(t)={\nu}_{\min, \mathrm{i}}(t)\cdotp {\mu}_{\mathrm{i}}(t)\cdotp \frac{L_{\mathrm{i}}(t)}{S_{\mathrm{la},\mathrm{i}}} $$

(26)

We note however that there is a large variability in the residual nitrogen in the senesced tissues, because cultivars can differ in the extent to which they retain chlorophyll after anthesis and, therefore, nitrogen in their leaves (Thomas and Howarth 2000). Nevertheless, replacing ν_{min, i}(t) with ν(t) in equation 26 (i.e., assuming no retranslocation) has almost no effect on the aboveground dry weight, leaf area index, and yield (data not shown).

Initialization of state variables

Each M³ simulation starts at sowing date t_s, defined with reference to the species of the intercrop that is sown the earliest. At this date, each of the state variables is set equal to 0, with the exception of the soil nitrogen amount at the sowing date N_s(t_s), which is set equal to an initial amount of nitrogen N_{s, 0}:

$$ {N}_{\mathrm{s}}\left({t}_{\mathrm{s}}\right)={N}_{\mathrm{s},0} $$

(27)

Once the crop emerges (S_i(t) = S_e,i, t = t_e.i), the other state variables are set to values larger than 0. The crop height is set equal to an initial value H_{0, i}:

$$ {H}_{\mathrm{i}}\left({t}_{\mathrm{e},\mathrm{i}}\right)={H}_{0,\mathrm{i}} $$

(28)

M³ assumes that each plant of crop species i has a certain leaf area per plant A_{0, i} (m² leaf plant⁻¹) at the time of emergence. The initial leaf area at the day of emergence is calculated from this leaf area and the plant density (plant m⁻² ground) ρ_i:

$$ {L}_{\mathrm{i}}\left({t}_{\mathrm{e},\mathrm{i}}\right)={A}_{0,\mathrm{i}}\cdot {\rho}_{\mathrm{i}} $$

(29)

M³ also assumes that the aboveground dry matter of crop species i at the day of emergence consists entirely of leaves. Therefore, the aboveground dry matter is calculated as:

$$ {B}_{\mathrm{i}}\left({t}_{\mathrm{e},\mathrm{i}}\right)=\frac{A_{0,\mathrm{i}}\cdot {\rho}_{\mathrm{i}}}{S_{\mathrm{la},\mathrm{i}}} $$

(30)

The initial amount of nitrogen in the crop as:

$$ {N}_{\mathrm{c},\mathrm{i}}\left({t}_{\mathrm{e},}\right)=\frac{A_{0,\mathrm{i}}\cdot {\rho}_{\mathrm{i}}}{S_{\mathrm{la},\mathrm{i}}}\cdot {\nu}_{\max, 0} $$

(31)

Field trial data

We used data from three experiments to calibrate and validate M³ for simulating growth and partitioning of light and nitrogen in mixtures of wheat. Experiment 1 and 2 were used to calibrate M³, while experiment 3 was used for validation. The experiments were carried out in different years in Wageningen, the Netherlands (51^° 59′ 24” N, 5^° 39′ 0″ E) and provided within-season observations of leaf area index and total aboveground dry matter. More details on these experiments can be found in Supplementary Tables S1 (Boons-Prins et al. 1993; Kropff 1989), S2 (Gou et al. 2017a; Gou et al. 2016), and S3 (Wang et al., unpublished data). From experiment 1, we used the control treatment monocrops of faba bean (cv. Monica) as monocrop (Boons-Prins et al. 1993; Kropff 1989). From experiment 2, we used the monocrop treatment for spring wheat (cv. Tybalt) (Gou et al. 2017a; Gou et al. 2016). Experiment 3 consisted of both monocrops and intercrops of spring wheat (cv. Nobless) and faba bean (cv. Fanfare) (Wang et al., unpublished data). Both species were grown within 1.5 m wide strip widths, with 6 rows per strip. During each experiment, irrigation was applied to avoid drought stress. The total amount of nitrogen applied as fertilizer was 20 kg N ha⁻¹ for faba bean in experiment 1. In experiment 2, 159 kg N ha⁻¹ in 2013 and 172 kg N ha⁻¹ in 2014 for wheat. In experiment 3, 125 kg N ha⁻¹ for wheat monocrops and 20 kg N ha⁻¹ for faba bean monocrops. In the intercrop treatment in experiment 3, the total amount of nitrogen that was applied to the wheat strips was 125 kg N ha⁻¹ land occupied by wheat strips and 20 kg N ha⁻¹ land occupied by the faba bean strips. Hence given that each crop covered half of the area, the overall amount of nitrogen applied to the intercrop is therefore (125 + 20)/2 = 72.5 kg N ha⁻¹ intercropped land.

Meteorological data

Minimum daily temperature, maximum daily temperature and daily global solar radiation data were obtained from the weather station De Veenkampen, west of Wageningen, maintained by the Wageningen University and Research chair group Meteorology and Air Quality. Supplementary Figs. S1 - S5 display these meteorological observations for the periods covering experiment 1 (1985–1986), 2 (2013–2014), and 3 (2019). Additionally, we used daily minimum and maximum temperature data from 1980 to 1982 from the same weather station to calculate the temperature sum between sowing and emergence for faba bean.

Model calibration and validation

Most crop parameters in M³ are adopted from literature (Table 1). The sub-sections below describe the estimation of the radiation use efficiencies (η_i), parameters related to crop height growth (H_0,i, r_h,i, H_max,i), and various phenological parameters (S_e,i, S_a,i, S_m,i, d_4,i) using data from experiments 1 and 2. We conclude the section by describing how M³ was validated using data from experiment 3.

Table 1 List of parameters, their symbols, values, units, and their meaning

Estimation of light use efficiency

A key parameter in M³ is the light use efficiency. This parameter was estimated by first approximating the time course of the leaf area index by linear interpolation between the measured values. We then used Beer’s law for each day between emergence and anthesis to calculate the daily amount of intercepted radiation, assuming an extinction coefficient of 0.6 (Zhang et al. 2014) for both wheat and faba bean. We then fitted a linear model to describe the relationship between the observed aboveground dry matter and the cumulated intercepted radiation. The radiation use efficiencies η_f and η_w were calculated as twice the slope of this relationship, because only half of the global radiation is photosynthetic active (i.e., f_p = 0.5). For wheat, we only used the data observed before the date of anthesis, so that the estimate of η_w is not affected by leaf senescence. The data set of faba bean did not contain enough observations of aboveground dry matter and leaf area indices before anthesis to do this. Therefore, we used observations throughout the growing season to estimate η_f.

Estimation of crop height growth parameters

We fitted the numerical solution of equation 12 to the height observations from experiment 1 for faba bean and experiment 2 for wheat in order to estimate H₀, H_m, and r_h. For experiment 1, only the data from 1985 contained height observations. For this nonlinear regression analysis, we used the Nelder-Mead algorithm (Nelder and Mead 1965) within the optim() function in the programming language R (R Core Team, 2018).

Phenological parameters, estimation of initial leaf area and temperature sensitive leaf senescence rate

Phenological observations and meteorological data were used to calculate the temperature sums from sowing to emergence (S_{e, i}), emergence to anthesis S_{a, i}, and anthesis to maturity S_{m, i} (equation 3) for faba bean (based on experiment 1) and wheat (based on experiment 2). Although the emergence dates of faba bean were known, the sowing dates in experiment 1 were unknown (Boons-Prins et al. 1993; Kropff 1989). We therefore assumed in the simulations of the field trials in experiment 1 that the fertilizer was applied at the emergence date. Therefore, we calculated the temperature sum from sowing to emergence for faba bean, S_{e, f} based on the sowing dates and observed emergence data from another data set, in which the same cultivar was grown in the period 1980–1982 in Wageningen (Boons-Prins et al. 1993; Grashoff 1990a, 1990b). The parameter d_{4, w} was calculated from phenological observations in experiment 2. This was not possible for d_{4, f}, because the only available phenological information in experiment 1 was the date of anthesis. We therefore adopted a value from the WOFOST crop growth model (Boons-Prins et al. 1993).

We fitted M³ on the leaf area index observations that were measured before anthesis to estimate the initial leaf area per plant A_{0, i}. Subsequently, we fitted M³, with the previously estimated value for A_{0, i}, on leaf area index observations that were measured after anthesis to determine the temperature-sensitive leaf area senescence rate s_sen. For these non-linear regression analyses, we used the L-BFGS-B algorithm (Byrd et al. 1995) within the optim() function of the R programming language (R Core Team, 2018).

Model validation

We simulated both the monocrops and the intercrops from the field trials in experiment 3 (Table S3) and compared observed leaf area indices and the total aboveground dry weights and storage organ dry weights to simulated values.

Quantification of model performance

We considered two indices to evaluate the correspondence between calculations with and the calibration data from experiment 1 and 2 and the validation data from experiment 3. These indices are the root mean squared error (RMSE) and the mean bias error (MBE), calculated as (Salo et al. 2016):

$$ \mathrm{MBE}=\frac{1}{n}\sum \limits_{\kappa =1}^n\left({y}_{\kappa }-{\hat{y}}_{\kappa}\right) $$

(30)

$$ \mathrm{RMSE}=\sqrt{\frac{1}{n}\cdotp \sum \limits_{\kappa =1}^n{\left({y}_{\kappa }-{\hat{y}}_{\kappa}\right)}^2} $$

(31)

where y_κ is the average of observed values of a focal variable y at measurement day κ, ${\hat{y}}_v$ is the simulated value of the same variable, n is the number of measurement days. MBE quantifies to what extent M³ overestimates (MBE < 0) or underestimates (MBE >0) the field observation. RMSE is the square root of the mean square deviations of the model from the data.

Sensitivity analyses

Global sensitivity analysis of nitrogen input level

M³ calculates the yield as the dry weight of the storage organs at the end of the growing season. We ran a global sensitivity analysis to investigate how the yield of monocrops and intercrops of wheat and faba bean would be affected by different nitrogen input levels. We used the nitrogen input level of experiment 3 as reference point and summarized different nitrogen input levels by means of the ration f_fert of the simulated amount of applied nitrogen to the amount of nitrogen that was actually applied in experiment 3.For instance, if f_fert = 0.5, then 62.5 kg N ha⁻¹ was applied to the simulated wheat monocrop (while the actual total amount applied to wheat in experiment 3 was 125 kg N ha), 10 kg N ha⁻¹ was applied to the simulated faba bean monocrop (while the actual amount was 20 kg N ha⁻¹), and the same amount to the wheat and fab bean strips in the intercrop leading to an average of 36.25 kg N ha⁻¹ for the intercrop (actual amount was 72.5 kg N ha⁻¹), f_fert was varied between 0 and 1.25 in steps of 0.025. We determined the final yields of wheat and faba bean in the monocrop and the intercrop for each value of f_fert. From these yields, we calculated the total yield of the wheat-faba bean intercrop, LER (equation 1) and TOI (equation 2).

Local sensitivity analysis of crop traits

We ran a local sensitivity analysis for each model parameter p, which represent crop traits, listed in Table 1 to see how small changes of size Δp in the value of this parameter affect the total yield of the intercrop. We calculated the relative sensitivity (also called elasticity) of the yields of wheat and faba bean in the intercrop and of the total yield of the intercrop to parameter p as (Van der Werf et al. 2007):

$$ {e}_{\mathcal{i}}(p)=\frac{p}{Y_{\mathrm{c}}(p)}\cdotp \frac{Y_{\mathrm{c}}\left(p+\Delta p\right)-{Y}_{\mathrm{c}}\left(p-\Delta p\right)}{2\cdotp \Delta p} $$

(32)

where Y_c(p) is the yield for a given value of parameter p and the subscript c indicates whether Y_c refers to the total yield of the wheat-faba bean intercrop, the yield of wheat in that intercrop or the yield of faba bean in that intercrop. ${e}_{\mathcal{i}}(p)$ is the relative sensitivity of ${Y}_{\mathcal{i}}$ to small changes in parameter p. In all cases, $\Delta p=\frac{1}{1000}\cdotp p$. We conducted this analysis for three levels of fertilization: the nitrogen input levels of experiment 3 (f_fert = 1), no fertilization (f_fert = 0) and an intermediate nitrogen input level (f_fert = 0.5).

Results

Model calibration using data from monocrops

When applied to monocrops, M³ underestimates the crop heights of wheat in both 2013 and 2014 during the early growth stages but was is in line with the observations in the later growing season (Fig. 5a). There was a good agreement between the measured and simulated values for crop height for faba bean (Fig. 5b). The model performed well in simulating the aboveground dry matter of wheat (Fig. 5c). It also performed well in simulating the aboveground dry matter of faba bean (Fig. 5d). There was good agreement between simulated and measured storage organ dry matter for both species (Fig. 5e and f), lending support to the use of M³ to simulate yields.

Fig. 5

Simulated versus measured values for crop height (a-b), total aboveground dry matter (c-d), and storage organs (e-f) for wheat (a,c,e) and faba bean (b,d,f). The solid lines are the one-to-one lines. RMSE and MBE are the root mean squared error and the mean bias error (eqs 30 & 31) . The data for faba bean were from experiment 1, conducted in 1985 and 1986. The data for wheat were from experiment 2, conducted in 2013 and 2014

Model validation using monocrop and intercrop data

We simulated both the monocrops and the intercrops from the field trials in experiment 3 (Table S3). Figure 6 compares measured and simulated results. There was a general good agreement between measurements and simulations, as indicated by the relatively low values for RMSE and MBE for both the monocrop (subscript “MC”) and the intercrop (“IC”). However, the model substantially overestimated the leaf area index for faba bean after anthesis (Fig. 6b). The effect of this overestimation on the total aboveground biomass of faba bean and the storage organ weights was limited, because also the measured leaf area indices were rather high after anthesis (over 4) indicating a near total interception of light. As a consequence, the modelled light interception and thereby the biomass was not much overestimated.

Fig. 6

Measured (dots) and simulated (lines) leaf area indices, aboveground dry weights and storage organ weights in experiment 3 (2019), grown in monocrop (MC) and intercrops (IC). RMSE and MBE are the root mean squared error and the mean bias respectively

Sensitivity analyses

Global sensitivity analysis nitrogen input level

We conducted a global sensitivity analysis of the yields of wheat and faba bean in monocrops and in intercrops to the nitrogen input level f_fert. Additionally, we have done this analysis for the total yield of the intercrop. Figure 7 shows the response of these yields to the nitrogen input level.. The yield of a monocrop of wheat was higher than the yield of faba bean or the total yield of an intercrop if f_fert ≥ 0.975. The simulated total intercrop yield was the higher than those of both monocrops (transgressive overyielding) if 0.775 < f_fert ≤ 0.975 (shaded area in Fig. 7) The simulated yield of a monocrop of faba bean was the highest if f_fert < 0.775.

Fig. 7

Global sensitivity analysis of various types of yields. The highlighted area represents the nitrogen input level for which transgressive overyielding occurs; i.e. where the transgressive overyielding index (TOI) is larger than 1

We also derived LERs (Fig. 8a) and the TOIs (Fig. 8b) from the simulated yields. TOI was larger than 1 if the nitrogen input level f_fert was between 0.775 and 0.975. This trend differed for the land equivalent ratio. It equaled 1.40, if f_fert = 0. It decreased with f_fert until it equaled 0.99 at f_fert = 1.25.

Fig. 8

Global sensitivity analysis for the land equivalent ratio and the transgressive overyielding index. The highlighted area represents the nitrogen input level for which transgressive overyielding occurs (i.e., TOI > 1)

Local sensitivity analysis

We conducted a local sensitivity analysis for three types of yields (wheat yield in intercrops, faba bean yield in intercrops, and total intercrop yield) at three different model nitrogen input levels (low: f_fert = 0, intermediate: f_fert = 0.5, high: f_fert = 1) to all model parameters. Figure 9 shows for all yield types to which 10 parameters these yields are the most sensitive at each nitrogen input level. At high nitrogen input levels, the wheat yield was most sensitive to the extinction coefficient of wheat k_w and radiation use efficiency of wheat η_w. At intermediate and low nitrogen input levels, the wheat yield was most sensitive to the maximum nitrogen concentration ν_{max, 0, w}. Nevertheless, k_w and η_w were still among the top 10 parameters to which the wheat yields were the most sensitive. The top 10 parameters to which the faba bean yield was most sensitive were identical for all nitrogen input levels. For all nitrogen input levels, the top 5 includes three phenological parameters for faba bean. These are the temperature sums from emergence to flowering (S_{a, f}) and from sowing to emergence (S_{e, f}), and the development stage above which no more biomass partitioning to other organs than storage organs and takes place (d_{4, f}). Additionally, this top 5 includes the radiation use efficiency of faba bean (η_f) and the maximum height of faba bean (H_{m, f}).

Fig. 9

Local sensitivity analyses of the yields of wheat (a-c), faba bean (d-f) grown in a wheat-faba bean strip intercrop and the yield of the whole intercrop (g-i) to model parameters under three different nitrogen input levels (f_fert ∈ (1,0.5,0)). For each figure, only the top 10 parameters for that yield (i.e for which the absolute value of the relative sensitivity is the highest) are shown. White bars indicate that increasing the parameter will increase the corresponding yield. Black bars indicate that increasing the parameter will decrease the corresponding yield. The last subscript of each parameters indicates whether it refers to wheat (‘w’) or faba bean (‘f’). Table 1 lists the meaning of the various parameter symbols and their symbols

Regardless the nitrogen input level, η_f, S_{a, f}, and d_{4, f} were always in the top 4 of parameters to which the total yield of the intercrop was most sensitive. η_w was also in the top 10 at high (rank 1) and intermediate input levels (rank 8).

Interestingly, none of the parameters related to height growth (initial height: H₀, maximum height: H_m, or height relative growth rate: r_h) of wheat or faba bean are in the top 10 most relevant parameters affecting the total intercrop yield under high and intermediate nitrogen input levels. Only at low nitrogen input levels, H_{m, w}, H_{m, f}, and r_{h, f} are present in the top 10. Increasing H_{m, w} will then decrease the total yield, while increasing the other two parameters will increase the total yield. In contrast, H_{m, w}, H_{m, f}, and r_{h, f} were always in the top 10 of parameters to which the faba bean yield was most sensitive. At high nitrogen input levels, these parameters ended in the top 10 parameters to which wheat yields is sensitive as well.

Discussion

Here we present a new minimalist mixture model for growth and resource partitioning in cereal legume intercropping, M³, and we calibrate the model to represent the growth and resource competition between wheat and faba bean in strip intercropping. The model performed well in simulating the aboveground dry matter, storage organ weights and the leaf area index of wheat, with the exception of substantially overestimating the leaf area index of faba bean after anthesis. One possible explanation is that M³ does not consider accelerated leaf senescence due to shading (Ackerly and Bazzaz 1995; Bazzaz and Harper 1977; Woledge 1972). Another explanation is that, although the model links senescence to temperature (equation 8), it does not explicitly consider the effect of short periods of high temperatures on senescence. Experiment 3 was conducted during the spring and summer of 2019. That summer had often higher maximum temperatures than average summers in Wageningen (Supplementary Fig. S5). The maximum daily temperature exceeded 30 °C on 11 days in 2019, while it exceeded 30 °C only two days in 1986 and it never exceeded 30 °C in 1985, i.e., the temperature was lower in the years on which the calibration data were collected than in the year 2019 in which the model was tested and which was characterized by one of the hottest summers on record (KNMI 2019). Finally, a further possible explanation is that the cultivar ‘Monica’, which was used for calibration (experiment 1), senesces later than anthesis more (Thomas and Howarth 2000) than the cultivar ‘Fanfare’, which was used for validation (experiment 3). This highlights the possibility of asymmetric cultivar influence on crops grown in pure or mixed stands, which can be explored by using models such as M³, especially when they are calibrated and validated with data from intercrops using same and different cultivars.

M³ assumes that the species in the intercrop take up nitrogen from the same pool of soil mineral nitrogen, similar to the assumption of another model that species in a strip intercrop take up water from the same pool (Tan et al. 2020). This simplifying assumption was made to keep the model simple in terms of the number of parameters and state variables and of the data requirements for calibration. In contrast, explicit modelling of these soil mineral nitrogen and root density gradients would require the addition of a 2-D model for root distribution and possibly growth, soil water movement, and mass flow and diffusion of soil mineral nitrogen (Chen et al. 2020). Also, the validation of such a model would require measurements of horizontal and lateral root distributions of both species (Chen et al. 2017; Gao et al. 2010; Li et al. 2017; Streit et al. 2019; Yu 2016), which are laborious to collect. Similarly to the YIELD-SAFE model (Van der Werf et al. 2007), M³ also makes simplifying assumptions that the rooting depth is constant throughout the growing season and that species in the intercrop have the same rooting depth. These various simplifying assumption were made to keep the number of parameters low. However, the current version of the model can therefore not be used to investigate to what extent differences in rooting depth and density affect interspecific competition and under which circumstances these differences can lead to niche complementarity for nutrient uptake. For instance, it has been reported that faba bean has a lower root density than wheat (Streit et al. 2019) and maize (Li et al. 2006), while faba bean has been reported to also have greater rooting depths than wheat (Bargaz et al. 2015). In order to test to what extent the assumptions regarding a constant rooting depth or one single pool of soil mineral nitrogen affects the distribution of nitrogen between two species in the intercrop, M³ should be validated on additional strip intercrop field trials with various fertilization levels and, if necessary, extended with a module to simulate the rooting depth, vertical and horizontal root growth, and root length.

M³ currently assumes that crops are grown under well-watered conditions, which is not always applicable. Although a model for strip intercrops grown under water-limited growth conditions has been published for wheat-maize strip intercrops (Tan et al. 2020), M³ does include the possibility to simulate nitrogen limited growth, which is of particular importance in cereal-legume intercrops. To the best of our knowledge, there is currently no model available that can simulate strip intercrops under both water and nitrogen limited conditions. Future contributions will aim at extending, calibrating, and validating M³ with a water balance and considering the effects of drought and water lodging on biomass production.

The effect of nitrogen input level on yields of wheat-faba bean intercrops.

We conducted a global sensitivity analysis for the nitrogen input level, f_fert. The results of this analysis show that the land equivalent ratio equals 1.00 for the highest simulated nitrogen input levels. In the unfertilized case, i.e., f_fert = 0, the land equivalent ratio, LER, was 1.4, indicating a major advantage of intercropping at this low level of nitrogen input (Fig. 8a). However, in absence of fertilization, the total yield of the intercrop was substantially lower than that of a monocrop of faba bean in the absence of fertilization. Several meta-analyses used LER to assess the performance of cereal-legume intercrops in terms of yield (Martin-Guay et al. 2018; Yu et al. 2015; Yu et al. 2016). A LER larger than 1 indicates that the yield of an intercrop of two species is higher than the yield of two monocrops of these species grown at the same location, in the same growing season. However, our results illustrate that this index alone does not provide information on whether the intercrop yield is higher than the yield of any of the species grown as monocrop (i.e., transgressive overyielding). This additional index TOI would provide additional insight in performance of the two systems. Figure 8 shows that transgressive overyielding only takes place within a relatively narrow intermediate range of nitrogen input levels and Fig. 7 provides an explanation for the relation between nitrogen input level and transgressive overyielding. At the highest fertilizer input level, the highest yield can be obtained by wheat monocrops due to the absence of nitrogen limitation. Consequently, there is neither substantial intraspecific nor interspecific competition for nitrogen and there is no mechanism for niche complementarity. At low nitrogen input, the highest yield is obtained by faba bean monocrops. Indeed, faba bean is able to fix N₂ from the atmosphere in order to fulfill most of its nitrogen demand, so the nitrogen input level has little effect on the yield, except when there is very little soil mineral nitrogen in the soil due to low mineralization rates. In contrast, the growth of wheat in the wheat-faba bean intercrop is substantially affected by nitrogen stress under low nitrogen input levels and the presence of faba bean cannot fully compensate the yield loss of the wheat. If specifically aiming at quantifying a yield advantage in legume-cereal intercropping systems, it is thus necessary to consider these separately for different nitrogen input levels and analyse explicitly whether transgressive overyielding is achieved. The latter can be achieved by including the transgressive overyielding index (Yu 2016) (equations 1–2), as we illustrated in Fig. 8b. Although the sensitivity analysis of the transgressive overyielding index shows for which nitrogen input levels the intercrop outperforms monocrops of both species in terms of yield, it should be mentioned that there are also indicators for other measures of crop performance. One example is the gross margin ratio, which quantifies whether the intercrop outperforms both monocrops in terms of gross margin (Van Oort et al. 2020). Another example is the fertilizer N equivalent ratio which indicates whether or not an intercrop that receives the same amount of fertilizer as two monocrops has a higher yield than two monocrops (Xu et al. 2020). The use of these indicators should depend on the aims of analyses and may result in different predictions of optimal nitrogen input levels in intercrops. Also crops are grown to support a complex food system in which multiple crop products are used and needed anyway. In that sense when wheat and faba bean need to be grown for a food system going for the most productive sole wheat would not produce the required faba beans and thus LER would be a relevant metric rather than TOI.

The effect of crop traits on intercrop yields and the interaction with nitrogen input level.

The relative sensitivities of the total wheat-faba bean intercrop yield to parameters related to height growth of any species (H₀, H_m, r_h) under high and intermediate nitrogen input levels were low in comparison to other parameters. In contrast, the individual yields of faba bean in the intercrop were substantially affected by H_{m, f}, H_{m, w}, and r_{h, f} under all nitrogen input levels. In addition, the yields of wheat were affected by these parameters, but only under high nitrogen input levels. This suggests that, under high and intermediate nitrogen input levels, increasing the height of one species will increase the light interception and accompanying biomass production of that species, but will also increase the extent to which this species shades the other species. As a consequence, the yield of this species decreases and the sum of the yields of both species will not be substantially affected. An explanation for the absence of an effect of wheat height on wheat yield at low nitrogen input levels is that, although wheat may intercept more light at these regimes, it can only convert this light to a limited extent to biomass due to nitrogen shortage. On the one hand, the total yield of the intercrop cannot be substantially affected by growing combinations of cultivars with different maximum plant heights under high and intermediate nitrogen input levels. On the other hand, the yield of one species can be favored, at the cost of a yield reduction of the other, by choosing cultivars of wheat and/or faba bean based on the relative heights that they can reach.

An interactive effect between species maximum height and planting geometry is also to be expected. Another model (Gou et al. 2017b) exploiting a similar concept of light interception applied to wheat-maize strip intercrops showed that the strip width can substantially affect the crop yields, as narrowing the strips increases the extent to which the tallest crop species shades the shorter crop species (Van Oort et al. 2020). Consequently, height-related parameters might have a more prominent effect on crop yields in intercrops with smaller strip widths and even more so in homogeneously mixed intercrops. This topic requires attention in future research, to evaluate the most suitable combinations of planting geometry and crop trait. M³ can be used for this purpose, even if the system of interest is an alternate row intercrop or a homogenously mixed intercrop. Specifically, M³ can simulate species in an alternate row intercrop by setting their strip widths equal to their row widths (Van Oort et al. 2020). And homogeneously mixed intercrop can be simulated as a system with strip widths equal to 0 (Gou et al. 2017a).

The local sensitivity analysis furthermore showed that some of the phenological parameters of faba bean (S_{f, f}, d_{3, f}) and light use efficiency of faba bean are important parameters for the total yield of the intercrop. At high nitrogen application rates, the light use efficiency of wheat also has a strong positive (e_{mc, w}(η_w) = 2.05) effect on wheat yield, while this effect is less important at intermediate nitrogen application rates (e_{mc, w}(η_w) = 0.97) or no nitrogen application (e_{mc, w}(η_w) = 0.2; not in the top 10 most relevant parameters). This result can be explained from the relationship between light use efficiency and leaf nitrogen concentration mentioned by Sinclair and Horie (1989) that reaches a plateau at high leaf nitrogen values where growth becomes purely light limited. While this concept was not directly built in M³, it is implicit in the incorporated concept of the critical nitrogen concentration (Greenwood et al. 1990). The crop biomass production is directly determined by light interception and light use efficiency at high crop nitrogen concentrations. Conversely, M³ includes a reduction factor of the light use efficiency (σ_{M, i}) at crop nitrogen concentrations below the critical concentration, which makes the light use efficiency of wheat a less dominant factor. It would be challenging to select good combinations of wheat and faba bean cultivars for increased light use efficiency, because light use efficiency is a complex trait, here assessed at a high level of aggregation. Indeed, light use efficiency is determined by various component traits at lower levels of aggregation. The component traits of light use efficiency include the photosynthetic capacity of leaves (Rodriguez et al. 1999), which is in turn determined by various component traits in the form of leaf biochemical parameters (Farquhar et al. 1980). In addition, the distribution of this photosynthetic capacity over the canopy is an important determinant of the light use efficiency (De Wit 1965; Goudriaan 1986; Kropff 1993; Rodriguez et al. 1999), and it is in turn affected by the distribution of nitrogen over the canopy of each species (Anten et al. 1995). If one would be interested in further improving the simulation of light use efficiency, models that are more mechanistic would have to be used. These models could elucidate which component traits of light use efficiency affect the yields in intercrops the most. However, this is only necessary for the components traits of the light use efficiency of wheat at high fertilizer application rates, as the light use efficiency of wheat had a less prominent effect on the yields at intermediate and low nitrogen input levels.

Conclusions

We developed a minimalist crop growth model for cereal-legume strip intercrops grown under nitrogen-limited, well-watered conditions. We calibrated and validated the model on monocrops and intercrops of wheat-faba bean. We aimed at identifying the key traits and nitrogen input levels leading to high yields of both species in the intercrop, high land equivalent ratios and the occurrence of transgressive overyielding. By means of global and local sensitivity analyses, we can conclude that:

• Transgressive overyielding only occurred in this system in a limited range of intermediate nitrogen input levels, whereas high values of land equivalent ratios occurred at low nitrogen input levels. Nevertheless, total yields of the crops and intercrops grown with low to intermediate nitrogen input levels did not exceed that of an optimally nitrogen-fertilized wheat monocrop.

• Height-related crop traits had little effect on the total yield of wheat-faba bean strip intercrops under high and intermediate nitrogen input levels.

• The maximum crop height of both species substantially affects the yield of faba bean and, at high nitrogen input levels, the yield of wheat in wheat-faba bean intercrops.

• The light use efficiency of faba bean can substantially affect the total yield of a wheat-faba bean intercrop, regardless the nitrogen input level. The light use efficiency of wheat has a substantial effect only at high nitrogen input levels.

• These results can provide guidance for choice and selection of cereal and legume varieties for sowing in intercrops. We recommend for further research to:

• Include transgressive overyielding indices in future meta-analyses that quantify yield advantages.

• Investigate the effect of crop height on the total yield of the intercrop in cereal-legume intercropping with narrow strip widths or no strips.

• Include mechanistic descriptions of radiation use efficiencies in crop growth models for intercrops, if one is interesting in increasing yields by changing their component traits.

• Further validate M³ by testing to what extent the assumptions of a single soil mineral nitrogen pool and a rooting depth that is constant over the season and equal for both intercropped species affect the nitrogen distribution between the two species in a strip intercrop.

• Extend M³ with a water balance and modules for water limited growth.