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

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. We developed a novel, parameter-sparse process-based crop growth model (Minimalist Mixture Model, M3) 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. 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. 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 1 and pLER 2 , 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: 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. Metaanalyses report average LER values for 1,32 for maizesoybean (Xu et al. 2020), 1.29 for intercrops with maize ), 1.17 for cereal-legume intercrops ) 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;: 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 3 − , NH 4 + ) from the soil to fulfil their nitrogen demands, legumes fulfill a large part of their nitrogen demand by fixing atmospheric nitrogen (N 2 ). 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 processbased 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. 2017Ritchie 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 CE-RES. In the case of APSIM, the model was found to have limited capabilities to simulate cereal-legume intercrops of wheat-faba bean 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 nitrogenlimited growth conditions. In this study, we developed, calibrated, and validated a process-based minimalist model (M 3 : 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 3 explicitly describes intra-and inter-specific interactions typical of strip intercrops and relies on a limited number of parameters, which have clear biophysical 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 3 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.

Model description
General overview M 3 M 3 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 2 leaf m /2 ground), ii) the total aboveground biomass (B i (t); kg dry matter m −2 ground), iii) the storage organ dry matter (Y i (t); kg dry storage organ matter m −2 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 −2 ground), and vii) the soil mineral nitrogen amount in the rooting zone (N s (t); kg mineral N m −2 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 3 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 3 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 3 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 3 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 3 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 3 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 3 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 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) 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 3 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.
Finally, M 3 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 3 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 3 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): 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: 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.

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: 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 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 2 leaf kg −1 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 −2 d −1 ), which is calculated as: 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 −1 , 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 0 (t) is the global radiation (J radiation m −2 land d −1 ) at time t. μ i (t) is the relative mortality rate of leaf biomass (d −1 ). 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 3 , 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: For the relative mortality rate μ i (t), M 3 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: Simulation of the aboveground biomass The rate of change of the aboveground dry matter is calculated as: which combines the production from equation 5 and the loss from equation 8. M 3 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: 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 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, Simulation of crop height M 3 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: where r h, i is the temperature sensitive relative growth rate of height growth (°C −1 d −1 ) 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 3 adopts a number of simplifying assumptions similar to the LINTUL3 model (Shibu et al. 2010). M 3 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 3 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: where f recov is the recovery fraction, i.e., the fraction of mineral nitrogen in the fertilizer that is not lost due to 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 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 −2 d −1 ). M is the mineralization rate (kg N m −2 d −1 ); we determined its value based on experimental data as the ratio of an assumed amount 25 kg N ha −1 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 −5 kg N m −2 d −1 . ∑ j U j t ð Þ is the sum of the nitrogen uptake rates of all species in the intercrop system (kg N m −2 d −1 ). The calculation of uptake rates will be specified in the next section. The change in crop nitrogen amount is simulated as: where 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 3 ( Fig. 1): they are all expressed in kg N kg −1 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 3 calculates the crop nitrogen concentration ν i (t) as: In order to calculate ν max, i (t), ν crit, i (t), and ν min, i (t), M 3 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: where are the β crit, i , β max, i , and β min, i are dimensionless shape exponents. ν max, 0, i (kg N kg −1 aboveground dry matter) is the maximum crop nitrogen concentration at low biomass (in this study, we assumed a biomass less than 1000 kg ha −1 ) 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 ℓ i þℓ j ℓ i equals 1. Figure 1 shows the nitrogen dilution curves schematically.
The nitrogen demand D i (t) (kg N m −2 d −1 ) of each species i in the system is calculated as: where t c is a relative rate (d −1 ) equal to 1 d −1 .
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 3 assumes that each species can fix a fraction f fix, i of its demand from N 2 in the atmosphere (Wolf 2012). M 3 calculates rate of nitrogen fixation as: 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: M 3 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 3 H 2 O m −2 soil) of species i is calculated as (Allen et al. 1998;Tan et al. 2020): 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 0 (t) is the reference evapotranspiration rate (m 3 H 2 O m −2 ground d −1 ) 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: 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: Each day the reduction factor of biomass production is calculated as a piecewise linear function of the plant nitrogen concentration at that time: 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 3 assumes that crops lose nitrogen via leaf senescence while the senescence of other plant organs is negligible. M 3 considers that the plants m i n i m i z e t h i s n it r o g e n l o s s b y n i t r o g e n retranslocation from senescing leaves to other plant organs, so that the nitrogen concentration of senesced leaves is ν min, i (t). As such, M 3 calculates the rate of nitrogen loss μ N (see equation 14) as: 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 3 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 : 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 : M 3 assumes that each plant of crop species i has a certain leaf area per plant A 0, i (m 2 leaf plant −1 ) 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 −2 ground) ρ i : M 3 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: The initial amount of nitrogen in the crop as: Field trial data We used data from three experiments to calibrate and validate M 3 for simulating growth and partitioning of light and nitrogen in mixtures of wheat. Experiment 1 and 2 were used to calibrate M 3 , 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  . 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 −1 for faba bean in experiment 1. In experiment 2, 159 kg N ha −1 in 2013 and 172 kg N ha −1 in 2014 for wheat. In experiment 3, 125 kg N ha −1 for wheat monocrops and 20 kg N ha −1 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 −1 land occupied by wheat strips and 20 kg N ha −1 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 −1 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 3 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 3 was validated using data from experiment 3.

Estimation of light use efficiency
A key parameter in M 3 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 0 , 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 w a s g r o w n i n t h e p e r i o d 1 9 8 0 -1 9 8 2 i n Wageningen (Boons-Prins et al. 1993;Grashoff 1990aGrashoff , 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 3 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 3 , 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): where y κ is the average of observed values of a focal variable y at measurement day κ, b y v is the simulated value of the same variable, n is the number of measurement days. MBE quantifies to what extent M 3 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 3 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 −1 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 −1 was applied to the simulated faba bean monocrop (while the actual amount was 20 kg N ha −1 ), and the same amount to the wheat and fab bean strips in the intercrop leading to an average of 36.25 kg N ha −1 for the intercrop (actual amount was 72.5 kg N ha −1 ), 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 wheatfaba 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): 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 i p ð Þ is the relative sensitivity of Y i to small changes in parameter p. In all cases, Δp ¼ 1 1000 Á 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).

Model calibration using data from monocrops
When applied to monocrops, M 3 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 3 to simulate yields. 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.

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.
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.

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 ).
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 0 , 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 3 , 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 3 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 3 , especially when they are calibrated and validated with data from intercrops using same and different cultivars.
M 3 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  (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;, which are laborious to collect. Similarly to the YIELD-SAFE model ( Van der Werf et al. 2007), M 3 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 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 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 3 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 3 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 3 does include the possibility to simulate nitrogen limited growth, which is of particular importance in cereallegume 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 3 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 wheatfaba 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;. 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 2 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 0 , 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 3 can be used for this purpose, even if the system of interest is an alternate row intercrop or a homogenously mixed intercrop. Specifically, M 3 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 3 , 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 3 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 cereallegume 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 wheatfaba 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 3 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 3 with a water balance and modules for water limited growth.