Title: Insights on nitrogen and phosphorus co-limitation in global croplands from theoretical and modelling fertilization experiments

Crossed fertilization additions are a common tool to assess nutrient interaction in a given ecosystem. Such fertilization experiments lead to the definition of nutrient interaction categories: e.g. simultaneous co-limitation, single resource response, etc. (Harpole et al., 2011). While these categories are commonly used in literature, what each category implies However, the implications of such categories in terms of formalism of nutrient interaction modeling are not remains unclear. To this end, we developed a theoretical analysis of nitrogen (N) and phosphorus (P) fertilization experiments based on the computation of ratios between plant demand and soil supply for each nutrient. The theoretical analysis is developed following two mathematical formalisms of interaction: Liebig's law of minimum (LM) and multiple limitation hypothesis (MH). As results of the theoretical framework, we defined what values of the limitation the values of the limitation the corresponding between most Harpole categories andby each nutrient when considered alone in the control experiment (i.e. without additional nutrient supply) are required to make an ecosystem in a given category. We clarified which categories are compatible with each interaction formalism assumed (LM or MH). In particular, Wwe showed that synergistic co-limitation could occur even using Liebig'sLM formalism under certain conditions, as a function of the amount of N and P added in fertilization experiments. e.g. if the ecosystem is N-limited in the control and if the amount of N added in the fertilization experiment is enough to switch the ecosystem into P-limitation. We then applied our framework with global maps of soil supply and plant demand for croplands to achieve their potential yield. This allowed us to estimate the global occurrence of each limitation category, for each of the possible interaction formalism. We found that a true co-limitation could affect a large proportion of the global crop area (e.g. ~42% for maize) if multiple limitation hypothesisMH is assumed.


Abstract:
Crossed fertilization additions are a common tool to assess nutrient interaction in a given ecosystem. Such fertilization experiments lead to the definition of nutrient interaction categories: e.g. simultaneous co-limitation, single resource response, etc. (Harpole et al., 2011). While these categories are commonly used in literature, what each category implies However, the implications of such categories in terms of formalism of nutrient interaction modeling are not remains unclear. To this end, we developed a theoretical analysis of nitrogen (N) and phosphorus (P) fertilization experiments based on the computation of ratios between plant demand and soil supply for each nutrient.
The theoretical analysis is developed following two mathematical formalisms of interaction: Liebig's law of minimum (LM) and multiple limitation hypothesis (MH). As results of the theoretical framework, we defined what values of the limitation the values of the limitation the corresponding between most Harpole categories andby each nutrient when considered alone in the control experiment (i.e. without additional nutrient supply) are required to make an ecosystem in a given category. We clarified which categories are compatible with each interaction formalism assumed (LM or MH). In particular, Wwe showed that synergistic co-limitation could occur even using Liebig'sLM formalism under certain conditions, as a function of the amount of N and P added in fertilization experiments. e.g. if the ecosystem is N-limited in the control and if the amount of N added in the fertilization experiment is enough to switch the ecosystem into P-limitation. We then applied our framework with global maps of soil supply and plant demand for croplands to achieve their potential yield. This allowed us to estimate the global occurrence of each limitation category, for each of the possible interaction formalism. We found that a true co-limitation could affect a large proportion of the global crop area (e.g. ~42% for maize) if multiple limitation hypothesisMH is assumed.

Introduction
In global assessments of crop ecosystem productivity limitations by nutrients, nitrogen (N) and phosphorus (P) are sometimes considered independently (Peñuelas et al., 2013); or they are considered together but without focusing on how the interaction modulates the limitation (Mueller et al., 2012). N and P cycles interact strongly with different processes that are key to this coupling (Achat et al., 2016). The most commonly studied interaction is related to the limitation of plant growth by nutrients: an increase in organ biomass (mainly composed of carbon, C) requires a given amount of both N and P, to respect stoichiometrical constraints. The interaction between carbonC and nutrients is usually represented by C:nutrient ratios for each organ. Plant growth is assumed to be limited when the demand for nutrients, estimated from C:nutrient ratios and C available for potential growth, is not satisfied by the supply of nutrient taken up by the plant. Due to incomplete knowledge about the mechanisms at the basis of the interaction and how these mechanisms are combined when integrating spatial scales (plant organ, individual, community, ecosystem) (Ågren et al., 2012;Davidson and Howarth, 2007;Sistla and Schimel, 2012;de Wit, 1992), the characterization of multiple element limitation remains an open scientific question. Two formalisms are generally used: Liebig's law of the minimum (LM) or the multiple limitation hypothesis (MH). In LM, plants are assumed to be limited by a single nutrient at a time, while in MH, it is assumed that plants adjust their growth patterns and thus they are co-limited by multiple nutrients simultaneously (Ågren et al., 2012). The MH formalism thus assumes that plants will mine the least available nutrient by using other resources. For instance, plants or groups of species growing in an ecosystem with a P-poor soil will invest C and N in the root system (and potentially to fungal mycorrhizae that form symbioses with plant roots (Ryan and 4 50 55 60 Graham, 2018)) to access more P (Davidson and Howarth, 2007). Both formalisms could be considered as macro-properties that reflect the same plant adjustments processes but, depending on the conditions, those adjustments may lead to an emerging behaviour that verifies one or the other formalism (Ågren et al., 2012). The further the supply of an essential nutrient deviates from a conceptual and theoretical optimum stoichiometry of plants, the more plants will follow the LM formalism (Ågren et al., 2012). LM is commonly assumed in many studies and is for instance used in most large-scale models dealing with multiple nutrient limitations (Barros et al., 2004;Goll et al., 2012;Mueller et al., 2012).
One way to assess the current nutrient limitation empirically is to provide single applications of +N, +P and +NP and to measure the increase in ecosystem productivity as compared to a control trial without any application. Such experiments are usually called fertilization experiments. By definition (Harpole et al., 2011), there is a true NP colimitation when the ecosystem is observed to respond to combined N and P addition only, or to both N and P when added separately. The different categories of nutrient limitation are summarized in Harpole et al. (2011) and in Table 1. Fertilization experiments are common in natural ecosystems and meta-analysis of these experiments have provided a global picture of nutrient limitation in natural ecosystems Elser et al., 2007;Harpole et al., 2011). Results from recent meta-analyses have shown that a true co-limitation is found in 28-42% of the studies Harpole et al., 2011 ones are more or less promoted by the interaction formalism assumed. To this end, we provided a theoretical framework of N and P fertilization experiments based on the computation of ratios between plant demand and soil supply for each of the two nutrients. The theoretical analysis is developed for two mathematical formalisms of interaction (LM or MH). This allowed us to define, for each formalism, the corresponding between Harpole categories and the values of the limitation by each nutrient when considered alone. Then, we analytically investigated how the choice of formalism modifies the NP limitation.
Finally, we applied our framework to the case of nutrient limitations in croplands. The justification of this choice is twofold: first, nutrient limitation is a key question in croplands at the global scale. For instance, (MacDonald et al., 2011) showed that 30% of cropland are characterized by negative soil P budget. Mueller et al. (2012) showed that 70% of the cropland where potential yield is not achieved could close their yield gap by increasing nutrient inputs. Second, For croplands, experiments with single and crossed N and P fertilizer applications of fertilizer (as defining "fertilization experiments") are not soas common in croplandsas those for natural ecosystems. This prevents us from having a global picture of N and P limitation based solely on observations, contrary to what was done in natural ecosystems (Elser et al. 2007, Harpole et al. 2011. Indeed, in cropland, fertilization experiments are usually characterized by one single addition for N (e.g. Di Paolo and Rinaldi, 2008;Salvagiotti et al., 2008) while for P, the same amount of fertilizer is applied each year for decades in so-called , On the other hand and it is usually difficult to retrieve the application level before the experiment (e.g. Deguchi et al., 2017;Restelatto et al., 2017), which prevents an accurate definition of the control in these cases.long-term field experiments with crops responding both to the annual supply of fertilizer and to the cumulative effect on soil P availabilityare common in croplands 6 100 105 110 115 (especially for P, e.g. Bai et al., 2013). This makes difficult to decipher the contribution of each nutrient. Moreover, in such experiments, many P treatments are tested and but for a given treatment,; the same amount of fertilizer is applied each year for decades, which makes the limitation in the long-term trial somehow non-representative to nutrient limitation happening in the surrounding fieldsprecludes analysis of current limitation.
Crossed fertilization additions are difficult to decipher from multi-nutrients and repeated fertilizer applications, as usually performed in croplands. To our knowledge, no meta-analysis of NP limitation in cropland field trials exists, which prevents us from having a global picture of N and P limitation based solely on observations. When a single application is the focus of a study, it is usually difficult to retrieve the application level before the experiment (e.g. Deguchi et al., 2017;Restelatto et al., 2017), which prevents an accurate definition of the control in these cases. Thus, our analysis, based on a theoretical analysis, is particularly adapted to investigate nutrient limitations in cropland. We applied our framework on global spatially explicit computations of soil supply and plant demand of N and P for croplands to achieve their potential yield, in order to assess the occurrence of co-limitation in croplands for each interaction formalism.we provided a theoretical framework of N and P fertilization experiments based on the computation of ratios between plant demand and soil supply for each of the two nutrients. The theoretical analysis is developed for two mathematical formalisms of interaction (LM or MH). This allowed us to define the corresponding between Harpole categories and the values of the limitation by each nutrient when considered alone. We analytically investigated how the choice of formalism modifies the NP limitation. Potential yield is here defined as the yield achieved without limitations of water and nutrients and without pest/diseases. In this work, Finally, we applied our framework on global spatially explicit computations of soil supply and plant demand of N and P for croplands to achieve their potential yield, in order to assess the occurrence of co-limitation in croplands for each interaction formalism.

Theoretical framework
Based on a framework commonly used in global studies (Goll et al., 2012;Kvakić et al., 2018), we defined the limitation of a nutrient considered alone as the ratio of its soil supply (S) and the demand by the plant to achieve its potential biomass (D): where SX and DX correspond to the supply and demand of the nutrient X, respectively (in kgX/ha/yr) with X is in {N,P}.
Crossed fertilization experiments are a common tool to assess nutrient limitation on a given site. They correspond to changes in nutrient supply in different combinations from the control (E1): addition of N alone (E2), P alone (E3) or N and P together (E4) (Fig. 1).
Based on the above equations defining the limitations of N and P (Eqs 1 and 2, respectively), theses changes in nutrient supply translate into limitations of each nutrient for each experiment E as follows: with AN and AP corresponding to the increase of N and P soil supply following addition of 9 155 160 165 N and P, respectively.
In the above framework, each nutrient is considered alone while the two nutrients interact. An ecosystem is thus defined by its NP limitation, called RNP in the following.
Two formalisms of interaction have been here considered to compute RNP from RN and RP: multiple limitation hypothesis (called MH in the following, Eq.7) or Liebig's law of the minimum (LM, Eq.8): where Ei is the experiment i. In MH, the limitations when the nutrients are considered independently (RN and RP) are multiplied to compute the NP limitation while in LM, the smallest one is selected.
We analytically investigated to which extent the choice of the formalism has an effect on the value of RNP for a given ( In fertilization experiments, nutrient limitation is assessed by looking at the change in productivity (pro) according to the addition of P alone (pro+P), N alone (pro+N) or N and P together (pro+NP). pro is here not expressed in absolute change but expressed relatively to the potential productivity (i.e. without any limitation). Harpole et al. (2011) defined different categories of limitation when considering the two nutrients in interaction. Each category is entirely defined by: i) the character null or non-null of pro+N and pro+P and ii) the relationship between pro+NP and (pro+N+pro+P) (column 3 of Table 1). Following Harpole et al. (2011), co-limitation exists when the increase in productivity following the addition of N and P together is strictly greater than the sum of increases in productivity when each nutrient is added alone (i.e. Δ pro +NP >Δ pro +N + Δ pro + P ). Any co-limitation is defined as a synergistic relationship. A given co-limitation is in addition considered as true if the responses to +N and +P are either both equal to 0 (i.e. Δ pro +N =0 and Δ pro +P =0 , simultaneous co-limitation, category A in Table 1) or both non-null (i.e. Δ pro +N ≠0 and Δ pro +P ≠0 , independent co-limitation, category B).
Here, we assumed that the change in productivity following the addition of +N, +P or +NP is equal to the change in RNP following the nutrient addition, i.e.: where Ei is the experiment i ( Fig. 1). This is a key assumption in our approach based on two simplifications described in details in the following. First, Tthrough these equations, we assumed that the productivity of a given experiment is proportional to RNP and that the slope of this relationship is equal to 1. In fact, a slope equal to 1 is not necessary to develop the theoretical analysis described in Text S1. As mentioned before, Harpole categories are defined through i) the character null or non-null of pro+N and pro+P and ii) the relationship between pro+NP and (pro+N+pro+P). These definitions are true even if the productivity of each experiment (and thus the different pro) is divided by the same slope. We keep here a slope equal to 1 for the sake of simplicity. Second, Eq.9-11 also imply that the relationship of proportionality between the productivity and RNP is true for all values of RNP, in the range [0-1]. In reality, the productivity vs. limitation relationship is very likely asymptotic (e.g. Bai et al., 2013). Here, we may approach this Text S1 for both MH and LM.
We showed in particular that to belong to the category "independent co-limitation" (category B in Table 1) with MH formalism, an ecosystem has to be characterized by both RN and RP in (]0,1)[ (a parenthesis instead of a square reversed bracket used in an interval means here that the corresponding endpoint is excluded from the interval; e.g. R in [0,1) means 0≤R<1). All other categories (A, C-G) require at least one ratio equal to 0 or 1: e.g.
Categories E,F,G are defined by Δ pro +NP =Δ pro +N +Δ pro +P and we showed that this requires at least one ratio equal to 1 with MH formalism (Text S2).
We showed that the formalism LM cannot represent true co-limitation, except in the very specific category A (i.e. RP=RN≠1). We found that synergistic co-limitation alone (categories C and D) can occur with LM but to be in these categories, the amount of N (if the control is N limited) or P (if the control is P limited) added in the fertilization experiments should be large enough to remove the initial limitation.
Conclusions of this analysis are summarized in Table 1 (columns 4 for MH and column 6 for LM).

Computation of spatially explicit RN and RP
We computed spatially explicit maps of RN and RP in croplands (0.5° latitude x 0.5°l ongitude) based on the computation of nutrient demand and soil supply. We then applied the above described theoretical framework on these RN and RP values to classify each grid-cell according to Harpole categories for the two interaction formalisms.
The computation of supply and demand maps used to estimate RN and RP are described below. To summarize, plant nutrient demand is based on nutrient harvest index data from the literature combined with spatially explicit distribution of crop potential yield (Ypot) (Mueller et al., 2012). The soil N supply has been approached by using a soil N budget taking into account fertilizer (mineral and organic), atmospheric deposition, biological fixation, and losses by ammoniac volatilization (Bouwman et al., 2011a), while the soil P supply is assessed by a potential root uptake model that accounts for soil P diffusion and soil P legacy effects. Such supply and demand estimates are representative to a growing season time-scale. More details about the computation of each variable can be found in Table S1.
Following Kvakić et al. (2018), demands for N and P to reach potential yields (DN and DP, in kg(N or P)/ha/yr) were derived from the combination of i) fixed parameters related to distribution of carbon (C) and nutrients between the different plant organs at maturity found in the literature and ii) spatially explicit potential yield (Ypot, in kgC/ha/yr): XHI corresponds to the nutrient harvest index (i.e. the ratio between the nutrient content of grain and the nutrient content of shoot, no unit), HI is the harvest index (i.e. the ratio between the carbon content of grain and the carbon content of shoot, no unit), RSR is the root/shoot ratio (no unit) and X%,grain and X%,root are nutrient concentrations (kg(N or P)/kgC) for grain and root, respectively. Kvakić et al. (2018) has shown that a XHI-based method provides similar demand estimates compared with others based on the nutrient concentration of all plant organs or QUEFTS (Sattari et al., 2014). The definition of the parameters used in Eq.12 (XHI, HI, RSR, X%,grain, X%,root) is based on nutrient and C biomass of different plant organs. These definitions as well as values used in the study are given in Table S2. Spatially constant values are here used. In particular, as the aim of our study is to assess nutrient limitation, we used organ concentrations derived from field experiments in stressed conditions in a multitude of climatic and socio-economic environments (Van Duivenbooden (1992) and Table S2).
Details about the Eq.12 are given in Text S3.
The potential yield (Ypot) is provided by Mueller et al. (2012) in tons of dry matter per hectare. In Mueller et al. (2012), the world grid-cells are divided into climate bins, defined by different combinations of growing degree days and amount of yearly precipitation; and within a climate bin, the potential yield characterizing this bin is defined as the area-weighted 95 th percentile of the grid-cell observed yields.
The supply of P (SP, in kgP/ha/yr) corresponds to the sum of a potential root uptake and a prescribed fraction (called ) of the inorganic content of total P fertilizer applied in the year considered Kvakić et al. (2018). The potential root uptake is determined by soil P availability and monthly root length density, following some assumptions about P diffusion in soil (Text S4). Soil P availability is derived from the current global distribution of soil P, as in Kvakić et al. (2018). The global distribution of soil P was determined by 15 295 300 305 310 315 combining information on farming practices, soil P dynamics, soil biogeochemical background, climate effect on soil P dynamics, etc., as well as the past variation of some of these drivers . Thus, we explicitly considered the soil P legacy effect, as it has been shown to be an important process (Ringeval et al., 2014;Sattari et al., 2012). Root characteristics parameters, root biomass at harvest (derived from Ypot, RSR and HI) and seasonality in root biomass ( sub-Saharan Africa), N applied in previous years is generally neglected due to the higher lability of N than P. We follow the same assumption as it is common in global modelling approaches Conant et al., 2013;Lassaletta et al., 2014).
Mineralization of soil organic N was also neglected as under steady state conditions it is expected to be compensated by N immobilisation in soil microbial biomass. N leaching was also neglected as in Bouwman et al. (2011), it is assumed that N leaching concerns only the surplus of annual soil agronomic N budget and occurs at the end of the growing season.
We recognize that the use of constant parameters at the global scale in the computation of supply and demand is too simple (Sadras, 2006). This is in particular true with respect to plant adjustments to nutrient limitations (Colomb et al., 2007) which are susceptible to modify nutrient organ concentrations. Also, some agronomic management such as cultivar diversity across World regions is susceptible to modify parameters, in particular HI. However, both plant adjustments (Franklin et al., 2012) and the effect of cultivar diversity on allocation (Folberth et al., 2016) are difficult to represent at the global scale.
Considering grid-cells independently in our uncertainty analysis (Text S5) made these parameters artificially vary in space.
Each term (SN, DN, SP, DP) is spatially explicit at half-degree resolution. An uncertainty related to each term has been considered (Text S5). Maize, rice and wheat are considered in this study (see the crop-dependent terms in Table S1) and the ratios computed are representative of the year 2000. Only grid-cells for which RP and RN could be computed are considered, which determines the crop area and the global crop production considered in our study (Table S3). In the Main Text, a specific focus is made on maize because it is the most widespread across latitudes. Caveats of our approach are discussed in Section 4.  Fig. S1). Large differences are also noticeable in regions with high limitations of both nutrients, such as the Western Russian Federation and Ukraine.

Occurrence of Harpole categories
We computed the occurrence of each Harpole category by using conditions onin terms of RN and RP , as described in Table 1. We checked that these occurrences are equal to the valuesoccurrences found when: numericalmodelling fertilization experiments are performed, RNP are computed for each experiment (Eq.7-8) and Eq.9-11 are then applied.
The Iincrease of N and P supply (AN and AP) in fertilization experiments are here equal to 30kgN/ha/yr and 5kgP/ha/yr, respectively and are spatially homogeneous for all cropland around the World. While our theoretical framework was initially developed for productivity (Section 3.2), we applied it here to cropland yield, which is consistent with the assumption of fixed harvest index as described in Section 3.1.
With the formalism MH, we found that true co-limitation occurs in 41.7±0.6 % of the global crop area for maize, via independent co-limitation (category B in Table 1). This category is found in the USA, South America, the Western Russian Federation and Ukraine (Fig. 4a). As showed theoretically, to belong to that category a crop has to be characterized by both RN and RP in (]0,1)[. In our simulations, these conditions occur for 42% of the maize crop area.
Synergistic co-limitation alone (categories C and D) occurs for 6.7±0.3 % of the global maize crop area and this is only explained by serial limitation N (category C, dark blue in Fig. 4a): no serial limitation P was found in our numerical application. This can be explained by the fact that RP (contrary to RN) is never null in our simulations because of the soil P legacy taken into account in our approach . This also prevents simultaneous co-limitation (A) from being found. The occurrence of colimitation at the global scale varies between crops (41.7±0.6 % for maize, 32.5±0.4 % for wheat and 18.7±0.8 % for rice, not shown). Except for few regions (e.g. India), grid-cells where the three crops are grown belong to the same limitation category for all crops (not shown): the difference in occurrence of co-limitation between crops is mainly explained by the crop-specific global distribution.
As theory shows (Text S1), the formalism LM cannot represent true co-limitation, except in the very specific category A (i.e. RP=RN≠1), which is never encountered in our study ( Fig. 4b  can occur in more than 15% of the global maize area with LM. However, this number is sensitive to the amount of N and P added in the fertilization experiments (called respectively AN and AP in Fig. 1). E.g. a cropland which is initially P-limited is classified in the category D if the amount of P added (AP) is sufficient to remove the P limitation (i.e. the cropland becomes N limited); otherwise, it belongs to the category F (Table 1).

Discussion
Previous studies estimating the occurrence of co-limitation in natural ecosystems were based on fertilization experiments performed around the World and provide some insights into the best way to represent the NP interaction (Elser et al., 2007;Harpole et al., 2011). Studies reviewing such experiments (characterized by a single application of nutrient) are not available for cropland, so we chose another strategy by computing the occurrence of co-limitation for each interaction formalism. Our work also clarifies the mathematical conditions in terms of supply/demand ratios required to place an ecosystem into a category of nutrient limitation, as defined by Harpole et al. (2011). In particular, we found that synergistic co-limitation can occur with Liebig's law of the minimum under certain conditions that are functions of the amount of N and P added in fertilization experiments, as already suggested by Ågren et al. (2012). While Liebig's law of the minimum is based on the limitation by a single nutrient at a time, it allows synergistic co-limitation to happen, which could be counter-intuitive. We found that, if multiple limitation hypothesis is the most appropriate way to represent nutrient interaction, co-limitation should occur for ~50% of the maize crop area (42% of true and synergistic co-limitation + 7% for synergistic alone co-limitation). The occurrence of true co-limitation in croplands would be of a similar magnitude to those reported for natural ecosystems (28% in The percentage of true co-limitations found here is higher than those reported for natural ecosystems in Harpole et al. (2011) , 42% in (28%) but similar to values found in Augusto et al. (2017)  The occurrences of the different limitation categories that we provided are a function of the spatial distribution of RN and RP, as posited by our theoretical framework. However, these maps are prone to uncertainty due to simplifications in our modelling approach. As mentioned in Section 3.1, some simplifications are related to the use of constant parameters at the global scale in the computation of supply and demand while plant adjustments and some farming practices are susceptible to modify them. Global changes are also very likely modifying yield and grain composition (e.g. Long et al., 2006;Müller et al., 2014) and this effect was not considered in our study which does not simulate temporal changes in nutrient limitation.
First, we recognize that the use of constant parameters at the global scale in the computation of supply and demand is too simple (Sadras, 2006) in particular with respect to plant adjustments to nutrient limitations (Colomb et al., 2007) which are susceptible to modify nutrient organ concentrations. However, the aim of our study is to assess nutrient limitation and thus, we used organ concentrations derived from field experiments in stressed conditions (Van Duivenbooden (1992) and Table S2). Global changes are also very likely modifying yield and grain composition (e.g. Long et al., 2006;Müller et al., 2014) and this effect was not considered in our study which does not simulate temporal changes in nutrient limitation. Besides, it is worth nothing that considering grid-cells independently in our uncertainty analysis (Text S5) made these parameters artificially vary in space.
Other simplification is related to the We use ofd potential yield provided by statistical methods based on maximum attainable yield within climate bins (Mueller et al. 2012) .
but sSuch approaches have difficulty distinguishing irrigated and rainfed crops and thus, the here used Ypot could be in fact water-limited in some places (van Ittersum et al., 2013b). Updates to the statistical methodology are ongoing to improve the separation between water-limited and irrigated yield potential (Mueller, personal comm.).
Alternative estimates of potential yield such as the ones simulated by Global Gridded Crop Models are very likely prone to huge uncertainties too (Müller et al. 2017). We also did not consider some agronomic managements that are susceptible to modulate nutrient limitation. In particular, we did not consider cultivar diversity across World regions. Such diversity is susceptible to modify parameters (in particular, HI) which are considered constant in space in our approach. However, cultivar diversity is difficult to consider at the global scale and up to now, it was mainly investigated through spatial variability in phenological development (van Bussel et al., 2015). To a much lesser extent, the effect of cultivar diversity on allocation (e.g. through variability in harvest index, susceptible to modify the here compute nutrient demand) was taken into account (Folberth et al., 2016). Also, some effects of crop rotation on nutrient limitation were not considered in our study. E.g. crop rotation can modulate the soil P availability because of difference in the strategies developed for enhancing nutrient acquisition among crops (Redel et al., 2007) and this effect was neglected here. However, N fixation by leguminous that can be incorporated within the crop rotation with cereals was indirectly considered in our study: our computation of N supply was not a function of crop (Table S1) and thus, the N supply budget encompasses an term for N fixation by leguminous occurring in the same grid-cell as cereals (Bouwman et al., 2011).
In our approach, the limitation of potential yield is computed by considering current farming practices to derive the supply. Current practices could be influenced by other limiting factors: e.g. if a crop is water limited, farmers can adapt their practices and reduce their nutrient applications accordingly. Sensitivity tests where the demand would be derived from actual yield (instead of potential yield, as in the Main Text) could help in the determination of areas where other limiting factors might play a role (Fig. S2). The next step is to consider more limiting factors together.
Our theoretical analysis has also few caveats. In particular, we assumed a linear relationship between RNP and the productivity of each experiment (Eq. 9-11). As 24 500 505 510 515 underlined in the method section, our conclusions are still valid if we assumed a linear relationship up to a value thresh if thresh replaces 1 in the definition categories given in Table 1. The value thresh is nevertheless theoretical because the calculated nutrient limitation (RN, RP, RNP) has no physical meaning and is disconnected from physical measure of e.g. soil P content (Olsen P, etc.). The fact that the transition between linear and plateau regimes occurs for the same RNP (1 or thresh) globally should be an acceptable assumption as we took into account the spatial variation in soil properties to compute the soil nutrient supply.
In our analysis, we computed Our computation of RN and RP, i.e. the increase in RN and RP required to increase RNP up to 0.75 and assessed how the choice of the interaction formalism has an effect on RN and RP. is based on the minimum "physiological" needs for plants. Behind the multiple limitation's mathematical formalism, an increase in RNP can be achieved for different combinations of increases in N and P (i.e. for different couples (RN, RP)): despite non-substitution at the molecular or cellular level (Sinclair & Park, 1993), one element can partly compensate for the other at the plant scale. Here, we considered only one couple (Fig. 3), while external variables such as the price or the ease of access to fertilizers will also influence the farmer's choice and could make him/her select another NP combination. This should be taken into account in future attempts to make link with scenarios of nutrient management and policy more straightforward. RN and RP could be translated to increase in soil supply by considering nutrient demand in each grid-cell. However such change in supply cannot be easily translated into a change in fertilizer, since our supply estimates take into account some processes occurring after the fertilizer application: for P, we take into account the dynamics of P in soil (diffusion and root uptake) while for N, we allow for NH 3 25 525 530 535 540 545 volatilization. Our nutrient requirement calculation is driven solely by nutrient limitation, independently of yield gap, contrary to previous estimates based on: soil quality indicators (with no distinction between N and P) (Fischer et al., 2012;Pradhan et al., 2015), statistical relationships between fertilizer application and yield (Mueller et al., 2012) or "N uptake gaps" based on yield gap and minimal/maximal values of the physiological N efficiency in aboveground biomass derived from the QUEFTS model (ten Berge et al., 2019;Schils et al., 2018). More generally, our nutrient limitation is not straight connected to the yield gap because the actual yield is not used in our computation. It is interesting to note that our computation of RN and RP is based on the minimum "physiological" needs for plants. Behind the multiple limitation's mathematical formalism, an increase in RNP can be achieved for different combinations of increases in N and P (i.e. for different couples (RN, RP)): despite non-substitution at the molecular or cellular level (Sinclair & Park, 1993), one element can partly compensate for the other at the plant scale. Here, we considered only one couple (Fig. 3), while external variables such as the price or the ease of access to fertilizers will also influence the farmer's choice and could make him/her select another NP combination. This should be taken into account in future attempts to make link with scenarios of nutrient management and policy more straightforward.
Our theoretical analysis has also few caveats. In particular, we assumed a linear relationship between RNP and the productivity of each experiment (Eq. 9-11). As underlined in the method section, our conclusions are still valid if we assumed a linear relationship up to a value thresh if thresh replaces 1 in the definition categories given in  (Ågren et al., 2012). It was also stipulated (Farrior et al., 2013) that  Tables   Table 1 (two pages). Nutrient limitation categories defined in Harpole et al. (2011) and occurrence for each crop in our modelling approach with MH formalism. pro+X is the change in productivity following the application of +X (with X=N, P or NP) in fertilization experiments. In the 1 st column, the y-axis defines ecosystem productivity and the dots correspond to the different experiments (white: control, blue: after addition of P, red: after addition of N, magenta: after addition of NP). Each category is defined as function of i) the character null or non-null of pro+N and pro+P and ii) the relationship between pro+NP and (pro+N+pro+P) (3 th column). Synergistic co-limitation means that Δ pro +NP >Δ pro +N + Δ pro + P . The different categories (columns 1-3) are derived from Harpole et al. (2011) while category B is restricted here to the "super-additive case" (subadditive or additive are neglected because they cannot happen in MH or LM, see Text S1). Long      A given grid-cell is defined by its (RN, RP) in the plan characterized by the base ( ⃗ R N , ⃗ R P ) . For a given gridcell and a given formalism, we called ⃗ u the shortest vector linking (RN, RP) and the curve (or segments) defining RNP=0.75. We called x and y the compounds of ⃗ u in the basis ( ⃗ R N , ⃗ R P ) , i.e. ⃗ u= ( x y ) . We defined Δ R N =max(0, x) and Δ R P =max (0, y) . In the above figure, two grid-cells are provided as an example: (RN=0.2; RP=0.5) for the black dot, and (RN=0.9; RP=0.1) for the black star. The formalism of interaction defines the (RN,RP) couples that make RNP=0.75: the blue curve defines RNP=0.75 for MH while the two orthogonal red segments define RNP=0.75 for LM. ⃗ u is provided for each grid-cell and each formalism (blue arrow for MH ; red arrow for LM). We explicitly plotted the RN and RP for the black dot and the two formalisms (solid black lines). Note that for the gridcell symbolized by the black star, Δ R N =0 for LM.  Table 1 and in Harpole et al. (2011) for MH (a) and LM (b) for maize. Category A corresponds to simultaneous colimitation, category B to independent co-limitation (super-additive), categories C and D to serial limitation (N and P, respectively), categories E and F to single-resource response (N and P, respectively) and category G to no response. For LM, whether one grid-cell belongs either to category C (dark blue) or to category E (cyan) depends on the value of AN. The same reasoning applies for categories D (yellow) and F (red) with AP.

Supporting Information
Supporting Text Text S1. Analytical characterization of the categories defined in Harpole et al. (2011) Text S2. Demonstration of ( R P =1 or R N =1)⇔(Δ pro +NP =Δ pro +N +Δ pro +P ) with the MH formalism Text S3. Computation of the nutrient demand (DN and DP) Text S4. Computation of the potential P uptake Text S5. Global values and uncertainty Text S6. Spatial distribution of RN, RP, RNP Text S7. Relationship between RNP and yield Text S8. Characterization of each category defined in Harpole et al. (2011) in terms of values for RP and RN with the MH formalism Supporting Tables  Table S1. Description and computation of the different terms used in Eq.1-2 of the Main Text Table S2. Parameters used to estimate the N and P demands (DN and DP, respectively) Table S3. Global crop area and production provided by global datasets and considered in our study Table S4. Values possible for RN and RP and the implications with MH formalism Table S5. For all crops, global values of supply (S), demand (D) and supply/demand ratio (R) for N and P when the two nutrients are considered as independent Figure S1. Spatial distribution of RN and RP Figure S2. The effect of using the real yield (instead of potential yield) on the computed nutrient limitation Figure S3. Grid-cell distribution in percentiles of different variables Figure S4. For maize, spatial distribution of RN and RP when N and P are considered as independent: average and standard-deviation of the 1000 replicates Figure S5. For maize, the spatial distribution of nutrient limitation when N and P are considered to be independent (bivariate plot of RN and RP) Figure S6. For maize, spatial distribution of RNP: average and standard-deviation for both formalisms of interaction Figure S7. Scatterplots of the ratio Yreal/Ypot provided by Mueller et al. (2012) vs. the simulated RNP at the country scale for maize