Nutrient Limitation Reduces Land Carbon Uptake in Simulations with a Model of Combined Carbon, Nitrogen and Phosphorus Cycling

Terrestrial carbon (C) cycle models applied for climate projections simulate a strong increase in net primary productivity (NPP) due to elevated atmospheric CO 2 concentration during the 21st century. These models usually neglect the limited availability of nitrogen (N) and phosphorus (P), nutrients that commonly limit plant growth and soil carbon turnover. To investigate how the projected C sequestration is altered when stoichiometric constraints on C cycling are considered, we incorporated a P cycle into the land surface model JSBACH (Jena Scheme for Biosphere–Atmosphere Coupling in Hamburg), which already includes representations of coupled C and N cycles. The model reveals a distinct geographic pattern of P and N limitation. Under the SRES (Special Report on Emissions Scenarios) A1B scenario, the accumulated land C uptake between 1860 and 2100 is 13 % (particularly at high latitudes) and 16 % (particularly at low latitudes) lower in simulations with N and P cycling, respectively, than in simulations without nutrient cycles. The combined effect of both nutrients reduces land C uptake by 25 % compared to simulations without N or P cycling. Nutrient limitation in general may be biased by the model simplicity, but the ranking of limitations is robust against the parameterization and the inflex-ibility of stoichiometry. After 2100, increased temperature and high CO 2 concentration cause a shift from N to P limitation at high latitudes, while nutrient limitation in the tropics declines. The increase in P limitation at high-latitudes is induced by a strong increase in NPP and the low P sorption capacity of soils, while a decline in tropical NPP due to high autotrophic respiration rates alleviates N and P limitations. The quantification of P limitation remains challenging. The poorly constrained processes of soil P sorption and biochemical mineralization are identified as the main uncertainties in the strength of P limitation. Even so, our findings indicate that global land C uptake in the 21st century is likely overestimated in models that neglect P and N limitations. In the long term, insufficient P availability might become an important constraint on C cycling at high latitudes. Accordingly, we argue that the P cycle must be included in global models used for C cycle projections.


Introduction
The neglect of nutrient limitation in carbon (C) cycle models has been criticized for being unjustified from an ecological point of view (Reich et al., 2006) and for overestimating terrestrial C sequestration in the future (Hungate et al., 2003; Published by Copernicus Publications on behalf of the European Geosciences Union.D. S. Goll et al.: Nutrient limitation reduces land carbon uptake Thornton et al., 2007).The incorporation of a nitrogen (N) cycle into the biogeochemical components of Earth System Models (ESM) generally reduces the CO 2 fertilization effect and C losses due to soil warming (Zaehle and Dalmonech, 2011).It even may change the sign of the feedback between climate and the terrestrial C cycle from positive to negative (Thornton et al., 2007;Sokolov et al., 2008).However, this is not generally the case (Zaehle et al., 2010a).In the current generation of coupled terrestrial carbon and nitrogen (CN) models, N limits productivity mainly in high latitude ecosystems (Sokolov et al., 2008;Thornton et al., 2007;Zaehle et al., 2010b), whereas the low latitude ecosystems, which are the most productive ones (Lieth, 1972;Field et al., 1998;Pan et al., 2011), are less N-limited.Low latitude systems are expected to show the highest direct CO 2 response of NPP (Hickler et al., 2008), and in climate change projections they are responsible for a substantial positive feedback between climate and the C cycle (Friedlingstein et al., 2006;Raddatz et al., 2007).At the same time, low latitude ecosystems are expected to be more P-limited (Townsend et al., 2011).Therefore, an inclusion of the terrestrial P cycle into global C cycle models seems essential to appropriately determine land C uptake.
Globally, in addition to N, P is the nutrient most commonly limiting plant growth and soil C cycling (Vitousek and Howarth, 1991;Aerts and Chapin, 2000).Theoretical (Walker and Syers, 1976) and observational (Reich and Oleksyn, 2004) studies suggest that at present, productivity of tropical forests and savannahs is generally constrained by P, while N constrains the productivity in temperate and boreal regions.The geographical patterns of N and P limitations are commonly explained by the age of the soils (Walker and Syers, 1976;Vitousek et al., 2010).Young soils tend to have a high P availability, but low N availability, because N is nearly absent in the parental material and accumulates over time from the atmosphere by biological N fixation.P enters the ecosystem by weathering of primary minerals.As soils age, the parental material can become P-depleted and P gets progressively occluded in secondary minerals.Therefore, old highly weathered soils tend to have low P availability.As glaciation has periodically reset soil development at high latitudes, these soils tend to be much younger than the ones at low latitudes, where glaciation was absent for hundreds of millions of years.
Additional mechanisms, operating on timescales of years to centuries (Vitousek et al., 2010), may cause the occurrence of nutrient limitation.Mechanisms limiting the availability of N include constraints on biological N fixation (Houlton et al., 2008) and losses by leaching and volatization, which may even occur to some extent in N deficient ecosystems (Davidson et al., 2007).As soil P is less mobile than N, losses are of minor importance for the occurrence of P limitation.A mechanism that affects both N and P limitation is the sequestration of N or P in an accumulating pool (Vitousek et al., 2010); in particular, the sequestration of nutrients in long-lived but non-photosynthetic pools of organic matter (Luo et al., 2004).Additionally, an enhanced supply of nutrients other than P, for example N (Perring et al., 2008), can cause P limitation.As the P cycle has a much slower turnover rate than the N cycle, such limitation can remain for centuries.Because mechanisms underlying N and P limitation differ from each other, the responses of N and P limitation to changes in atmospheric CO 2 , climate, and land-use change likely differ, too (Davidson et al., 2007;Vitousek et al., 2010;Gruber and Galloway, 2008).
Despite evidence of widespread phosphorus limitation on plant productivity from fertilizer experiments (Elser et al., 2007) and modelling (Wang et al., 2010;Zhang et al., 2011), to our knowledge, so far no study has assessed to what extent land C uptake projected by C-only models can be realized when stoichiometric constraints of N and P on the build-up of organic matter are accounted for.Results from Zhang et al. (2011) suggest that during the 20th century the global C balance can be fully quantified with just N, while at continental scale the P cycle is also of importance for the C balance.However, the uncertainties in simulating P cycling are large (Wang et al., 2010;Zhang et al., 2011).So far, a detailed analysis of whether and how P limitation may develop in the future has not been done.
As it is not yet clear where and how strong N and P constrain today's plant productivity and soil C turnover (Elser et al., 2007;Zaehle and Dalmonech, 2011;Townsend et al., 2011), we developed a new modelling concept of nutrient limitation to avoid prescribing nutrient limitation in the inital model state.We assume that during thousands of years of stable Holocene climate and relatively low atmospheric CO 2 concentration, plants have adapted to their environment such that their growth is limited by multiple resources (Field et al., 1992;Townsend et al., 2011).We thus hypothesize that for present day, ecosystems are co-limited by the availability of N and P, which is in broad terms consistent with a meta-analysis of N and P manipulation experiments across global biomes (Elser et al., 2007).As nutrients are needed for the build-up of organic compounds (C sink) from carbohydrates, nutrient uptake and C fixation (C source) have to be adjusted such that C sinks and sources are in balance.The adaptation of photosynthesis to nutrient availability is to some part represented in photosynthesis models, as the parameters of these models are specific for each plant functional type (PFT).Based on global data compilations, Kattge et al. (2009) concluded that the maximal rate of carboxylation (V c max ) is affected by N and P availability.Reich et al. (2009) showed that leaf P, an indicator of P availability, modulates the N dependence of photosynthesis.The nutrient limitation hidden in the parameterization of photosynthesis is called background nutrient limitation (BNL).
Today's atmospheric CO 2 concentration is unprecedented in the last 2 million years (Luethi et al., 2008;Hoenisch et al., 2009), and it is projected to increase further with a rate unseen in the geological past.In C-only models, photosynthesis is independent from the actual need of carbohydrates, and therefore the projected increase in plant productivity is invalid if the carbohydrate sink can not grow as fast as photosynthesis increases.To account for such a possible imbalance, plant productivity must be adjusted to the actual need of carbohydrates, which can be derived from the stoichiometry of plant tissue.The difference between the potentially possible productivity under elevated CO 2 and the sink corrected productivity is defined as the effect of CO 2 -induced nutrient limitation (CNL).CNL is per definition absent in the preindustrial state.In aggregate, the modeling approach taken herein does not attempt to quantify the overall background nutrient limitation, and instead estimates how much less C can be stored on land when changes in N and P supply are simulated; our approach is distinct from one which would quantify total N and P limitation as is done with fertilizer experiments (Elser et al., 2007) and in terrestrial ecosystem models (Sokolov et al., 2008;Thornton et al., 2007;Zaehle et al., 2010b;Wang et al., 2010;Zhang et al., 2011).
In this paper we aim to quantify how much of the land C uptake projected by a C model can be sustained when accounting for stoichiometric constraints on C storage in organic matter.Further, we aim to identify the processes responsible for the occurrence of P limitation.To do so, we implemented a P cycle into the biogeochemical cycle model (Parida, 2010) of JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg) (Raddatz et al., 2007), which includes C and N cycling.The model was applied under the SRES (Special Report on Emissions Scenarios) A1B scenario with and without N and/or P cycles to allow a separation of the effects of N and P limitation on land C cycling.In addition, a set of sensitivity experiments was performed to identify the main drivers of nutrient limitation.

Model
JSBACH describes the processes for land hydrology, phenology, and biogeochemical cycling (Fig. 1).The submodel for the biogeochemical cycles of JSBACH is driven by net primary productivity (NPP), leaf area index (LAI) and climate variables (soil temperature, soil moisture, runoff) obtained from full climate simulations with the Max Planck Institute Earth System Model (MPI-ESM) (Roeckner et al., 2011).By driving the submodel for biogeochemical cycling with NPP and LAI, we omit effects of CO 2 -induced changes in N or P availability on these properties because, so far, it is unclear whether we can generalize about how elevated CO 2 affects the interactions between nutrients and leaf level photosynthesis (Reich et al., 2006).
Background nutrient limitation (BNL) is indirectly considered in the model as part of the PFT-specific parameterization of the C cycle, which is based on measurements in present day ecosystems and therefore reflects the mean of present day nutrient conditions for a single PFT (Knorr, 2000).However, the range of different soil environments on which a single PFT can grow may lead to a variability of BNL within a PFT, which is not captured by our approach.The additional nutrient limitation caused by the increase in atmospheric CO 2 (CNL) is computed dynamically in the model assuming globally uniform constant C : N : P stoichiometry for the different biosphere compartments, such as active and woody biomass and soil organic matter.The model further employs a supply-demand-approach in which plant nutrient demand is derived from the ratio of potential C-allocation to the plant's C : nutrient ratios, while nutrient demand by soil microorganisms is derived from the potential immobilization fluxes.Subsequently, the total nutrient demand by plants and soil microorganisms is compared against the supply.If the supply is insufficient to meet the demand, the demand is diminished.These strong simplifications allow us to simulate the N and P cycles with a minimum number of parameters.The merit of our approach is that it avoids the introduction of additional parameters which are hard to quantify, unlike more mechanistically oriented CN models (Xu-Ri and Prentice, 2008;Jain et al., 2009;Zaehle and Friend, 2010;Esser et al., 2011), and thereby reduces its dependence on model parameters which are difficult to constrain by available ecosystem observations.The assumption of stoichiometric inflexibility is a strong simplification and may lead to an overestimation of CNL, but we argue it is an appropriate first-order approach given the limited process understanding and data to quantify the full nutrient dynamics including adaptation at plant and ecosystem level (Reich et al., 2006).The implications of these assumptions on nutrient limitation are discussed in detail in Parida (2010).

The C cycle
The C cycle in the biogeochemical cycling part of JSBACH is driven by NPP derived from the photosynthesis scheme, following Farquhar et al. (1980) for C 3 andCollatz et al. (1992) for C 4 plants and including an explicit dependence of productivity on atmospheric CO 2 concentration (Knorr, 2000).NPP is allocated in fixed fractions to root exudates and the three plant pools: a wood pool representing lignified plant tissue, a active pool for active plant tissue (leaves, fine roots, etc.), and a reserve pool containing C stored as sugars and starch.Upon leaf shedding, C is transferred from the active pool and the reserve pool to the non-lignified litter pool.The wood pool has a fixed turnover time, and C is transferred to a woody litter pool which differs from the non-lignified litter pool in its turnover time.Microbial activity mediates the transfer of C from the litter pools to the atmosphere and to the slow soil pool.The latter pool represents slowly decomposing organic matter and is depleted by heterotrophic respiration.The respiration rates depend on soil moisture, nutrient availability, and exponentially on soil temperature.The soil compartment consists of one pool (slow) representing slow-decomposing organic matter.All carbon pools except the reserve pool have a corresponding nutrient pool with variable C : N : P ratio (slow, non-lignified litter) or fixed C : N : P ratio (rest).There is one mobile plant pool representing mobile nutrients stored internally in plants.Soil mineral nitrogen is represented by a single pool (soil mineral pool), while mineral P is represented by labile (available) pool and sorbed pool.The opposing triangles indicated that the flux is controlled by phosphorus (red triangles), nitrogen (blue triangles) or both (mixed triangles).

The N cycle
The N cycle is described in detail in Parida (2010).N is taken up by vegetation from a single soil mineral N pool and allocated to the woody and active plant tissue according to the N : C ratios of the these pools.Prior to leaf shedding a certain fraction of N is retranslocated to the mobile N pool, representing internal N storage.Upon leaf shedding, the remaining N is transferred to the non-lignified litter pool.The wood pool has a fixed turnover time, and N is transferred to a woody litter pool.Microbial activity mediates the transfer of N from the litter pools to the slow soil pool.As the N : C ratio of the slow soil pool is larger than the N : C ratios of the litter pools, additional N from the soil mineral N pool is transferred (immobilization).The slow soil pool represents slowly decomposing organic matter and is depleted by biological mineralization.Biological N fixation, the flux of N from the atmosphere to the soil mineral N pool, is a function of NPP, comparable to the approach by Thornton et al. (2007).In this case, we use potential NPP to account for the positive effects of increases in N demand and carbohydrate supply on N fixation (Vitousek et al., 2002).We account for the effect of P limitation on N fixation (Reich et al., 2006) by scaling the fixation rate using the P limitation factor (Eq. 2).Denitrification, the loss flux to the atmosphere, is simulated as a function of the soil mineral N pool (Meixner and Bales, 2003).N losses by leaching are described analogously to P leaching using the approach of Meixner and Bales (2003).

The P model
The P cycle in JSBACH is designed to diagnose the P distribution on land and to adjust C fluxes for P availability.To include the cycling of P, the same methodology as for N is employed (Parida, 2010).Assuming that the C allocation is not affected by nutrients, potential C allocation and heterotrophic respiration are computed.We refer to these fluxes as potential C fluxes.In combination with the stoichiometry of the involved pools, we use the potential C fluxes to derive diagnostic P fluxes: where F C i j is the potential C flux from pool i to pool j , F P i j the corresponding diagnostic P flux, and r i and r j the P : C ratios of the pools i and j , respectively.P is released if r i > r j or sequestered if r i < r j .The latter case may only be realized to the extent that additional P is supplied.Therefore, we introduced the following rule: If the diagnosed demand is larger than the supply, P fluxes are diminished according to the available P supply.The supply is given by the labile P (P l ) and the demand is given by the sum of plant demand (D veg ) and microbial demand (D micr ).Based on the assumption that vegetation and soil organisms are equally strong in their competition for P, we use a single limitation factor to diminish plant and microbial demand.The amount of potential NPP which cannot be allocated due to P limitation is interpreted as an increase in autotrophic respiration.This reflects Table 1.Model variables (Fig. 1).For the C cycle, the indexes i and j are: A (atmosphere), a (active), w (wood), r (reserve), la (non-lignified litter), lw (litter wood), and s (slow).For the P cycle, the indexes i and j are: a (active), w (wood), m (mobile), la (non-lignified litter), lw (litter wood), s (slow), l (labile), and sr (sorbed).The index e stands for the element (N or P).

Symbol Description units
the high cost of nutrient foraging in a nutrient limited system (Hogberg et al., 2003).This is supported by a large-scale fertilization experiment in a N-limited boreal forest (Olsson et al., 2005), where fertilization reduced autotrophic respiration while aboveground productivity increased.We define the P limitation factor (f P limit ) as where the term in square bracket is the maximum rate at which the labile P pool can supply P. Note that in the discretized formulation the labile pool can at most be depleted during a single model time step ( t).We thus set this maximum rate to P l t .If the N cycle is active in addition to the P cycle, the limitation factor of the more limiting nutrient is taken to correct potential C fluxes and their corresponding diagnostic nutrient fluxes.
where f N limit is calculated analogously to f P limit (Parida, 2010).This is the only modification in the structure of the C cycle in comparison to Parida et al. (2011).
The P cycle consists of nine pools (P i ; see Fig. 1 and Table 1).P enters the terrestrial ecosystems by the weathering of P-bearing minerals.Weathered P enters the labile P pool (P l ), which is available to plants and microbes.Labile P is rapidly adsorbed onto soil particles.The associated pool for sorbed P is denoted by P sr .There are three pools representing P in vegetation.One pool is for P in active tissue (P a ), one pool is for P in woody tissue (P w ), and one pool represents mobile P (P m ), which acts as an internal buffer of plant P. Plants and microbes take up P from the labile P pool (P l ).Upon litterfall, P enters the litter pool for woody litter (P lw ) and the pool for non-lignified litter and fastdecomposing (life time ≈1.5 yr) soil organic matter (P la ).P in slow-decomposing organic matter (life time ≈100 yr) is represented by one single pool (P s ).Three P pools have constant P : C ratios (r i ), namely the active pool (i = a), wood pool (i = w), and woody litter pool (i = lw).P pools with a constant r i can be derived from the corresponding C pools (C i ) by The other P pools have a constant r i for the incoming flux but not for the pool itself (P la , P s ), or have no corresponding C pool and thus no r i (P l , P m , P sr ).The mobile P pool (P m ) is a short term P storage internal to plants.It is filled by P retranslocated prior to litterfall and is depleted by P used to allocate as much of the potential NPP (NPP pot i ) to C w and C a as possible.
The fraction of NPP pot i , which is directly allocated by use of P m , is named NPP dir i .NPP dir w and NPP dir a are given by Litterfall and wood shedding connect plant pools with the litter pools.Non-lignified litter is made of fast decomposing organic matter and leaves.The dynamics of the non-lignified litter P pool (P la ) are where the first term describes the P influx from litter fall, the second term describes the flux of P from herbivore excrements which is not directly available to biota, and the third term arises from the P released by biological mineralization of litter.We assume that active plant material (P a ) is consumed by herbivores at a constant rate (ε) and immediately excreted (Parida, 2010).We separate the excrement into labile (β l ) and fast decomposing (1 − β l ) P compounds which enter the non-lignified litter pool (P la flux F C la s is a demand flux (r la < r s ), therefore it is multiplied with f limit .α i indicates the fraction of decomposed non-lignified litter which enters the atmosphere by respiration.The P of the respired fraction of F C la s is mineralized and enters the P la .
The change in the P content of slow-decomposing soil organic matter (P s ) is given by the net of immobilization and mineralization fluxes: The first two terms describe P immobilization due to decomposition of non-lignified litter and woody litter.Mineralization occurs as biological mineralization (d s P s ) and phosphatase-mediated biochemical mineralization (F BC ).
The decomposition rate d s is computed by a Q 10 model as described in the Appendix B. Biochemical mineralization is modeled using a rate constant M min if P limitation is absent (f P limit = 1), but increases under P limitation until a maximum mineralization rate (M max ) with a lag time of 30 days.There is no observational evidence to support 30 days, but the results of this study are insensitive to the time lag.If P limitation occurs, we use a 30 day mean of f P limit (f limit ) to increase the rate until the maximal mineralization rate is reached.
This formulation is based on the concept of McGill and Cole (1981) that biochemical mineralization is P demand driven.The biochemical mineralization flux, which is mediated by N-rich enzymes (Phosphatases), depends on the availability of N. If the N cycle is active in addition to the P cycle, we use the 30 day mean of the difference between f P limit and f N limit (f limit ) instead of f limit in Eq. ( 9).
Thereby we account for a decreased stimulation of biochemical mineralization in cases when N is limiting in addition to P. If N is limiting more strongly than P, a stimulation of F BC is inhibited.If N is limiting, but less than P, the stimulation is less pronounced.The equations governing the dynamics of labile and sorbed P (P l and P sr ) are adopted from Wang et al. (2010).It is assumed that the labile P pool and the sorbed P pool are in equilibrium on a daily time step.The relationship between the amount of labile P and sorbed P is described by the Langmuir equation (Barrow, 1978) in its differential form: where S max is the maximum amount of sorbed phosphorus in the soil, and K s is an empirical constant.
Under these assumptions the dynamics of the labile P pool can be described by (see Appendix A for derivation) The inputs to the labile P pool are biological P mineralization (F M ), biochemical P mineralization (F BC ), dust deposition (F D ), P weathering (F W ), and directly available P in excrement of herbivores (εβ l P a ).Losses are leaching, strong sorption, and demand fluxes (f limit (D veg + D micr )).The loss rate of sorbed P to strongly sorbed and occluded P forms is assumed to be proportional to the amount of sorbed P in the soil ( 1 τ sr P sr ) (Wang et al., 2010).The flux of leached P is a function of runoff (R) and the amount of P l dissolved in solution (γ P l ).Biological P mineralization (F M ) is the sum of the mineralization flux from the slow pool and P leached from freshly sheded wood:

Nutrient demands
The P demand of vegetation (D veg ) is given by the amount of labile P needed for the allocation of the potential NPP to the vegetation compartments of wood (NPP The microbial P demand is the sum of all potential P immobilization fluxes, which are the fluxes from litter pools (non-lignified litter and woody litter) to the slow pool.When litter is decomposed, a certain fraction of its C is respired to the atmosphere, and the remaining C becomes slowly decomposing organic matter.Denoting the respired fraction by α i , the microbial P demand becomes: By respiration some P becomes available, which is subsequently immobilized by microbes, and therefore is subtracted from the microbial demand.The non-lignified litter pool does not have a fixed r i , instead the term P la C la is used in the above equations.The decomposition fluxes are and where i can be lw or la.d i is the decomposition rate of litter pool i, and its calculation is described in the Appendix B.

Parameter estimates
The parameters of the P model are globally uniform, except for the biochemical mineralization rates which depend on the soil order, and their values are shown in Table 2.We use information from the global data set of plant traits, TRY (Kattge et al., 2011), to compute the N : P ratios for the active and wood pools, as well as the for the ratio used for the computation of the flux from the active pool to the non-lignified litter pool.From these N : P ratios, P : C ratios are derived using the N : C ratios from Parida (2010).The P : C ratio of the woody litter is set to a value within the range of spread of P : C of bark (Kattge et al., 2011), as data for heartwood is scarce.The initial P : C ratio for the slow decomposing organic matter is taken from Smil (2000).
Leaching is calibrated so that, first of all, the global leaching flux compares well with estimates from empirical modelling (Seitzinger et al., 2005), and in addition, according to our concept of limitation, that CNL by P is absent under preindustrial conditions.
The biochemical mineralization flux is calibrated such that the N : P ratio of soil organic matter is close to the estimates by Yang and Post (2011).The outcomes of this study are insensitive to parameter values for the adjustment intensity (η) (not shown) of the biochemical mineralization flux in Eq. ( 9), which is set arbitrarily due to a lack of data.The dynamics of the labile and sorbed P are parameterized based on measurements on US soils using the USDA soil classification (Wang et al., 2010) because there is no more comprehensive data.Therefore, the history and the characteristics (soil texture, redox potential, etc.) of the different soil types reflect global patterns to the extent that US soil types are representative for the global diversity of soils.

Experimental setup
The simulation setup is adopted from the CN simulations performed by Parida (2010).We run the submodel of JS-BACH for the biogeochemical cycles independently from the rest of JSBACH driven by net primary productivity (NPP), leaf area index, soil temperature, soil moisture and surface runoff.The forcing variables are extracted from simulations performed with the MPI-ESM (Roeckner et al., 2011) under the framework of the ENSEMBLES project (Hewitt and Griggs, 2004).The MPI-ESM consists of the atmosphere model ECHAM5 (Roeckner et al., 2003) and its land surface scheme JSBACH (Raddatz et al., 2007) coupled to the ocean model MPI-OM (Marsland et al., 2003).The simulations by Roeckner et al. (2011)  19 vertical atmospheric layers, driven by prescribed atmospheric CO 2 concentration (SRES A1B scenario 1860-2100) and anthropogenic land-cover changes.Accordingly, our JS-BACH stand-alone simulations are conducted at T31 resolution, but without land-cover changes, using a map of potential vegetation instead.Output variables of the MPI-ESM are provided for every PFT at every grid independent of the prescribed vegetation cover.Therefore, the initial distribution of PFTs in simulations with the stand-alone model can be different from the one used in the forcing simulation.To investigate the long-term evolution of nutrient limitation, we performed an additional simulation with the full MPI-ESM to derive forcing data where we prolonged the original ESM simulation for an additional 400 yr keeping the atmospheric CO 2 concentration fixed at the level of 2100.In this period, climate slowly approaches a new equilibrium.
Using the submodel of JSBACH for the biogeochemical cycles, we brought the element cycles into equilibrium (less than 1 % change in the 30-yr mean of global pools) using a repeated 30-yr cycle of the forcing data from the coupled simulation (1860-1889).We equilibrated the element cycle consecutively in three steps.Firstly, we ran the model in the C-only mode reaching equilibrium after 5000 yr.Secondly, we prolonged this simulation by switching on the N cycle for an additional 5000 yr.Thirdly, the P cycle was brought into equilibrium after prolonging the CN equilibrium simulation for an additional 8000 yr in the CNP mode (submodel run with N and P).We made sure that this sequential spin-up procedure did not result in a different equilibrium state than a simultaneous spin up.During the spin up, N and P deposition were kept constant at the level of 1860 and present day, respectively, as the human impact on P deposition is marginal (Mahowald et al., 2008).
We conducted several transient simulations with the submodel: without nutrients (C-only), with N (CN), with P (CP), and with N and P (CNP) (Table 3).Land use, landuse change, and fertilizer application are not accounted for.N deposition is read in from maps (Galloway et al., 2004;Dentener et al., 2006) which are interpolated in time proportionally to the increases of atmospheric CO 2 concentrations of the forcing simulation until 2030, and afterwards extrapolated also using CO 2 of the forcing simulation, following Wang and Houlton (2009).Present day P deposition maps (Mahowald et al., 2008) are used for the whole simulation period.
To test the robustness of our results, we conducted several sensitivity simulations (Table 3).A set of simulations was performed with PFT-specific C : N : P ratios from Wang et al. (2010) and McGroddy et al. (2004) (CN pft , CP pft ).As a redistribution of P from unavailable pools to available ones can affect nutrient limitation, we varied parameters controlling the size of redistributable P, namely sorbed P (CP mo , CP ox ), and modified the adaptive strength of biochemical mineralization (CP bc ).In another simulation (CN dp ), we kept N deposition constant at the present level to test for the effect of increasing N deposition rates.Furthermore, we performed simulations with flexible stoichiometry, as well as a set of simulations with perturbed stoichiometric parameters.The latter simulations are described in detail in the Appendix C.

Present day: comparison with observations
For present day, the model simulates a terrestrial C stock of 2881 Gt(C), 24 % of it stored in vegetation (682 Gt(C)), 8 % in litter and fast decomposing soil organic matter (229 Gt(C)), and 69 % in slow-decomposing soil organic matter (1970 Gt(C)).When N and P cycling are taken into account (CN, CP, CNP) the present day C stocks are only slightly reduced (< 1 %).The simulated vegetation C storage is slightly above the higher end of inventory based estimates that range from 560 to 652 Gt(C) (Saugier and Roy, 2001).This relatively high vegetation C stock can be attributed to the land cover prescribed in our simulation, which differs from the real world as crop land is replaced with potential vegetation, which stores more C than crops.The simulated soil C stock is in the range of observation-based estimates that range from 1500 to 2300 Gt(C) (Post et al., 1982;Batjes, 1996), depending on the considered soil depth.
The decadal means of the net land to atmosphere fluxes for the 1980s, 1990s, and 2000s (Table 4) are between 2.2 and 3.7 Gt(C) yr −1 when nutrients are not considered.When N and P are considered, the fluxes are reduced up to 10 %.The simulated fluxes are higher but of the same magnitude as the ones by Friedlingstein et al. (2010).Note that we do not account for land use in our simulation, in contrast to Friedlingstein et al. (2010).When land use and land-cover change are explicitly accounted for, the land C uptake is reduced by 20 % to 30 % (Parida, 2010).The inclusion of the P cycle into JSBACH only slightly changes present day C cycling because CNL by P is mostly absent between 1860 and 2000 and reduces present day NPP by less than 3 % (not shown).Overall, the model performs well in simulating C Table 4. Simulated decadal land C uptake (Gt (C) yr −1 ) with and without accouting for nutrients compared to published estimates of the natural terrestrial sink from the Global Carbon Project (GCP) (Friedlingstein et al., 2010).Note that we omit land use (change) in our simulations.cycling given the large uncertainties in quantifying present day C cycling (Friedlingstein et al., 2010).For a detailed discussion of the C and N cycle, see Parida (2010).
The current state of the terrestrial P cycle is even less constrained by observations than the C and N cycle.Estimates of total terrestrial soil P range from 40 Gt(P) (Smil, 2000) to 200 Gt(P) (Jahnke, 1992).P stocks are often derived from C or N stocks using assumptions about the stoichiometry (Smil, 2000).Soil P is mostly measured in agricultural soils (Johnson et al., 2003), which are heavily influenced by fertilizer and manure application (Smil, 2000).Measurements of the fraction of soil P which is available to biota are prone to high uncertainties (Johnson et al., 2003).The mean sizes of the simulated P pools (CNP) for present day are summarized in Table 5, and the spatial distribution is shown in Fig. 2. When compared to modelling and empirical estimates, the simulated P storages are at the lower end of estimates.The upper ends of the estimates for P in soils of 200 Gt(P) (Jahnke, 1992) and in vegetation of 3 Gt(P) are criticized for being unrealistic (Smil, 2000;Wang et al., 2010).Our study supports these findings unless we have overestimated C : N ratios.However, the C : N ratios used in JSBACH are relatively low (Parida, 2010).The N : P ratios for the vegetation compartments are derived from a comprehensive stoichiometric data set for plants (Kattge et al., 2011), which is an improvement to earlier studies based on more restricted data sets.
The model seems to underestimate the amount of labile P in highly weathered soils when compared to a recent data compilation (Table 6) (Yang and Post, 2011).The compilation by Yang and Post (2011) includes the data which were used to calibrate the sorption model (Wang et al., 2010).Two aspects make the estimation of plant available P difficult.First, the interpretation of P measurements is not straightforward, as the measured fractions cannot be linked directly to physiological ones like plant available P (Johnson et al., 2003;Yang and Post, 2011).Second, the area based estimates by Yang and Post (2011) depend on assumptions about the bulk density of soils, which are highly variable within one class of soil.However, in the next section, we will show that the amount of labile P in the initial model state is of relatively small importance for the occurrence of P limitation.The annual P uptake by vegetation is shown in Fig. 2  and Table 7. From measurements, the annual P uptake can be estimated as the sum of P required to replace P in litterfall, P in wood increment, and P required for root growth, assuming equilibrium.Johnson et al. (2003) estimated the mean annual P requirement of temperate forests to be 0.52 ± 0.38 g(P) m −2 a −1 .For temperate broadleaf forests we simulate a mean (± standard deviation) uptake rate of 0.52 ± 0.23 g(P) m −2 a −1 .For tropical evergreen forests, P lost by litterfall, which makes up the largest part (approximately 60 %) of the P uptake in temperate forests (Yanai, 1992;Johnson et al., 2003;Yang and Post, 2011), is in the range of 0.2 to 0.4 g(P) m −2 a −1 (Vitousek, 1984;Valdespino et al., 2009;Yang and Post, 2011).Yang and Post (2011) report a higher mean loss rate of 0.6 g(P) m −2 a −1 using a litterfall database.The simulated mean P uptake of evergreen tropical forests is 0.68 ± 0.32 g(P) m −2 a −1 .Thus, despite spatially uniform CNP ratios, the model is able to simulate reasonable P uptake rates for forests.However, it may underestimate P uptake in grasslands.In the budget of two semiarid grasslands given by Cole et al. (1977), plant uptake was 0.53 and 1.2 g(P) m −2 a −1 , respectively.We simulate lower uptake rates of 0.24 ± 0.22 and 0.31 ± 0.28 g(P) m −2 a −1 for C 3 and C 4 grasslands, respectively.The low uptake rates are mainly caused by a low NPP of 340 and 420 g(C) m −2 a −1 for C 3 and C 4 grasslands.Observation-based estimates of grassland NPP range from 200-2000 g(C) m −2 a −1 , but may be even higher (Scurlock et al., 2002).For the 20th century, we can rule out an influence of low P availability on P uptakes rates via P limitation, as present day NPP is reduced by P limitation by less than 3 % in the model.In the sensitivity simulation with PFT-specific C : P ratios, P uptake by vegetation is reduced in tropical PFTs by 10-20 %, but increased in temperate PFTs by 40-50 %.Overall, the mean uptake rates do not improve using a PFT-specific set of C : P ratios.As the uptake is strongly controlled by the P demand of vegetation, the differences in P uptake rates between PFTs are dominated by the differences in NPP.
Immobilization of P by microbes is generally larger than the uptake by plants (Yang and Post, 2011).In a budget of two grassland sites, the immobilization fluxes exceeded the uptake by vegetation by factors of 2.6 and 4.7 (Cole et al., Table 6.Mean (± standard deviation) for simulated labile P (g(P) m −2 ) that typically supports forests compared to measurements (Yang and Post, 2011).The stocks, which are averaged over a 30-yr period , are from simulations with present day soils (CP), uniform low sorption soils (CP mo ), and uniform high sorption soils (CP ox ).1977).In our simulation, the immobilization fluxes of grasslands are globally 2.8 and 2.5 times larger than plant uptake for C 3 and C 4 grasslands, respectively.In the case of N, the immobilization flux is 1.7 times larger than plant uptake.The more dominant role of microbes in P immobilization compared to N is supported by the finding of Aponte et al. ( 2010) that 8.8 % of total P is immobilized by microbes compared to 4.7 % in the case of N.
The mineralization of P from slowly decomposing soil organic matter can occur by biological and biochemical mineralization (McGill and Cole, 1981).Biochemical mineralization is mediated by extracellular enzymes which specifically cleave out P from organic matter; thus it is controlled by the P demand of biota (Stewart and Tiessen, 1987).In contrast, biological mineralization is driven by the energy demand of microorganisms, as it is coupled to heterotrophic respiration.Although phosphatase activity, which is a qual-itative measure for biochemical mineralization, is common in soils (Stewart and Tiessen, 1987), the quantification of the mineralization rates in the field is not yet possible.The simulated biochemical mineralization fluxes are in the same order of magnitude as the biological mineralization fluxes (Table 7).In a simulation without biochemical mineralization, the amount of P in soil organic matter was 16.2 Gt(P)), which is 3.13 times more than in the standard simulations.The increased amount of P in soil organic matter in the simulation without biochemical mineralization is comparable in magnitude with the results of Wang et al. (2010) about the effect of biochemical mineralization on soil P.
Biochemical mineralization partly decouples P mineralization from decomposition of soil organic matter, namely N mineralization.Therefore, the N : P ratio of the slowly decomposing soil organic matter can vary in the model, and we used it to calibrate the strength of biochemical mineraliza- tion.The lower end of the N : P ratio is given by a prescribed C : N : P ratio, while the higher end depends on the strength of biochemical mineralization.The simulated N : P ratios of slow-decomposing soil organic matter of 20.3, 26.8, and 37.7 for slightly, intermediate, and highly weathered soils, respectively, compare well with data compiled by Yang and Post (2011).The agreement of simulated and observed N : P ratios indicates that in the model the biochemical mineralization is in a reasonable order of magnitude, unless the C : P ratio of newly formed organic matter used in the model is wrong.

Land C sink until 2100
In the C-only simulation, land accumulates 550 Gt(C) from 1860 to 2100.The uptake is in the range of land C uptake simulated by different C cycle models (Friedlingstein et al., 2006).As in Friedlingstein et al. (2006), we do not account for land-use change, which would reduce the land C uptake.The C sequestered on land is stored in vegetation (+315 Gt(C)), soil (+195 Gt(C)) and litter and fastdecomposing soil organic matter (+40 Gt(C)) (Fig. 3).When considering nutrient limitation, the projected land C uptake for the period 1860-2100 is reduced by 13 %, 16 % or 25 % due to limitation by nitrogen (CN), by phosphorus (CP) or by both (CNP), respectively (Fig. 3c).The strengths of nutrient limitations in the given model setup are not affected by the assumption of stoichiometric inflexibility (see Appendix C).
Parameter perturbation simulations reveal that the strength of N and P limitation depends to some extent on the parameterization of stoichiometry, but the ranking of limitations is not affected (Appendix C).P limitation reduces the amount of C stored in vegetation more strongly than N (Fig. 3d).This is explained by the geographic occurrence of N and P limitation (Fig. 4).P limits C uptake mainly at low latitudes, where the C taken up is predominately stored in vegetation (see Fig. D1 in Appendix D).At high latitudes the soil sink dominates.
The partly additive nature of N and P limitation on land C uptake can be attributed to their distinct geographic occurrence (Fig. 4).In addition, in the model, N limits mainly grasslands, while P limitation is nearly absent in grasslands (not shown).Hence, N and P limitation are to a large degree additive, despite their geographical co-occurrence in some regions.The lack of synergetic effects between the nutrient cycles can be attributed to the simple representation of N-P interactions in our model and may not reflect the actual effects.For example, N fixation is represented by an empirical function which cannot account for all pathways of N-P interactions underlying N fixation.To rule out that the pattern of nutrient limitations is a result of the globally uniform C : N : P ratios, we performed additional sensitivity simulations with a PFT-specific set of C : N : P ratios from Wang et al. (2010).In these simulations the pattern of limitations is not changed, although the strength of N and P limitation is slightly reduced (see Fig. D2 in Appendix D).
To test the influence of soil P sorption capacities on P limitation, we performed simulations with globally uniform parameters for the dynamics of labile and adsorbed P. The parameter for the sorption capacity (S max ) constrains the amount of adsorbed P in the initial state.As sorbed P is a buffer for labile P, the extent to which increased P demands can be sustained by a redistribution of adsorbed P to labile P is constrained by S max .The typical tropical soils, Oxisols and Ultisols, have the highest sorption capacities and the fertile soils of temperate regions, Alfisols and Molisols, have lower capacities.Therefore, as expected, an artificial exchange of all soils with Molisols (CP mo ) does not alleviate but increases P limitation drastically during the 21st century (Fig. 4).CNL by P becomes even larger than CNL by N at high latitudes.The low sorption capacities of soils in central Europe and eastern North America make these regions vulnerable to CO 2 -induced P limitation.An exchange with Oxisols reduces P limitation in our simulations due to the high sorption capacity.The accumulated land C uptake between 1860 and 2100 in CP mo and CP ox simulations is reduced by 19 % and 4 %, respectively, compared to the C-only simulation.When compared to the simulation CP, the land C uptake is less limited by P in the CP ox , despite the stocks of labile P being smaller (Table 6).This shows that the stocks of labile P do not necessarily correlate with P limitation.It is difficult to assess how realistic such a redistribution of adsorbed P to labile forms is.Evidence from short-term CO 2 enrichment studies suggest that a redistribution of P from unavailable to available forms does occur under elevated CO 2 (Johnson et al., 2004;Khan et al., 2008).However, it is unclear which processes are responsible for the observed redistribution and how these results from short-term experiments translate to natural ecosystems affected by a steady increase in CO 2 over centuries.
Biochemical mineralization is an additional mechanism that can shift P from unavailable to available forms.We assume that biochemical mineralization is highly adaptive to P limitation.In the CP simulation, biochemical mineralization rates increase by 37 % by the end of the 21st century.This is higher than an increase in biochemical mineralization by 18 % caused by doubling of CO 2 in the simulations by Wang et al. (2007).As the biochemical mineralization rates so far cannot be quantified from measurements, we conducted a sensitivity simulation (CP bc ) in which we kept the biochemical mineralization rates constant at the pre-industrial level.This results in a reduction of land C uptake by CNL by P of 39 % for the period 1860-2100, while the geographic pattern is only affected marginally (not shown).This sensitivity test shows, in addition to the sensitivity simulations CP ox and CP mo , the importance of adaptive processes which shift nutrients from unavailable forms to available forms for projecting nutrient limitation in a transient state.
High inputs by N deposition can shift ecosystems from N-limited to P-limited by increasing the P demand due to N-stimulated productivity and by acidifying soils which increase P sorption.In JSBACH, the effect of increasing N deposition rates on the geographic pattern of N and P limitation is marginal, as the differences in the latitudinal means of CNL between the simulation with increasing N deposition and the one with constant N deposition (CN dp ) are small (Fig. 4).In general, the effect of N deposition on the land C uptake during the 21st century is small in JSBACH owing to our concept of CNL (Parida, 2010).In addition, the P sorption characteristics are prescribed in the model.Hence, we may underestimate the effects of N deposition.While the simulations show high uncertainties in the quantification of nutrient limitation, a robust outcome from the simulations is the occurrence of CNL by P in low latitude and high latitude ecosystems during the 21st century and the nearly additive nature of P and N limitation.The geographical pattern of CNL by P and N is comparable to the pattern of total nutrient limitation in simulations by Wang et al. (2010).This pattern is to a large extent robust against the assumptions about biochemical mineralization and soil sorption, and it is not a result of special calibration procedures; this pattern automatically emerges in the simulation when atmospheric CO 2 concentration increases.

Land C sink after 2100
The global mean annual temperature over land increases in our simulations by more than 5 K during the 21st century (Fig. 3a), which is comparable to other models (Friedlingstein et al., 2006).Temperature is a main driver of the biogeochemical cycles.Under ongoing climate change, the cycling of elements gradually changes, leading to an imbalance between the C and the nutrient cycles.After stabilization of temperature and climate at the new level, biogeochemical cycles are slowly approaching a new equilibrium.To analyze nutrient limitation in the new climate state, we prolonged the simulation for 400 yr keeping CO 2 concentration constant at the level of 2100, but letting climate evolve.
Warming has two counteractive effects on the land C cycling in JSBACH.On the one hand, C losses by autotrophic and heterotrophic respiration increase.Concepts for their representation in models are currently under investigation.On the other hand, temperatures at high latitudes become more suitable for plant growth and the growing seasons extend, and therefore NPP increases.Overall, land becomes a source of C in JSBACH as the losses outweigh the sinks (Fig. 3c).Increasing temperature also directly influences the N and P cycling by enhancing mineralization in many parts of the world.The increase in mineralization of N and P by 35 % globally at the end of 2100 is comparable to the enhancement of N mineralization by 45 % measured in a 7-yr soil warming experiment with a comparable increase in temperature (+5 K) (Melillo et al., 2011).However, as temperature effects on N and P mineralization will strongly depend on interactions between temperature and soil water, such comparisons must be made with considerable caution.A decrease in tropical NPP at mean land temperatures warmer by as much as 7 K compared to today due to climate change is a feature of several coupled climate-carbon cycle simulations, although the amplitude of the drop in NPP varies considerably among models (Raddatz et al., 2007).A negative effect of increasing temperature on tropical productivity is supported by observational evidence (Clark, 2004).An increase in the NPP of temperate and boreal systems due to warming is projected in all simulations analyzed by Friedlingstein et al. (2006) (T. Raddatz, personal communication, 2011).
Part of this increase can be attributed to an advancement of leaf unfolding.In JSBACH, leaf unfolding of temperate broadleaf forests appears 17 days earlier at end of the 21st century than present day.The reported changes in late 20th century tree phenology correspond to an advancement of leaf unfolding of 16-24 days per K of warming, assuming that the average warming over the last five decades was about 0.6 K (Solomon et al., 2007).This relationship may not apply to the future, as a decline in freezing events may delay the breaking of dormancy, and therefore the effect of warming may be much smaller and could be even reversed (Morin et al., 2009).The phenology model of JSBACH accounts for this inhibitory effect, and the more moderate advancement of leaf unfolding is closer to the simulations by Morin et al. (2009), which project an advancement by 5 to 10 days during the 21st century under the SRES scenarios A2 (+3.2 K) and B2 (+1 K).
As a result of the increase in NPP, the N and P demands increase at high latitudes.In the tropics and subtropics, however, nutrient demands decrease as NPP declines with increasing temperatures.In JSBACH the increases in autotrophic respiration may be overestimated, as the model does not account for temperature acclimatization of autotrophic respiration (Wythers et al., 2005;Tjoelker et al., 2009) nor adjust for complexities of the instantaneous temperature response curve of respiration (Tjoelker et al., 2001).Therefore, the nutrient demand may be higher than our simulations suggest.In addition, we can not rule out that we may have underestimated tropical N limitation due to the simplistic representation of N fixation in the model, which does not account for possible negative effects of warming on N fixation in the tropics (Wang and Houlton, 2010).The stocks of soil organic matter decrease at low latitudes releasing nutrients, while at high latitudes the stocks of soil organic matter increase and nutrients get progressively sequestered (see Fig. D1 in the Appendix D).Together, the changes in demand and sequestration of nutrients result in an increase of CNL at high latitudes, and a decline in CNL at low latitudes (Fig. 5).As a consequence of the decline in tropical CNL, the land C stocks of these regions are nearly the same for all simulations (C, CN, CP, CNP) at the end of our simulations.
At high latitudes, N limitation declines while P limitation stagnates (Fig. 5).This shift from N to P limitation is more pronounced in the simulation with PFT-specific C : N : P ratios (CNP pft ) (see Fig. D2 in Appendix D), as P demands are higher for temperate and boreal forests due to lower C : P ratios.However, it is not known whether these lower C : P ratios merely reflect luxury consumption at present (Agren, 2008); if so, the P demand is overestimated in the simulation CP pft .Changes in species composition to more P efficient plants, which are not accounted for in our simulations, could also alleviate the increase in P demand.In general, the mechanisms underlying stoichiometric flexibility on all levels (leaf, single plant, ecosystem) are not yet fully understood to account for such flexibility in a comprehensive manner.In the model, the predominance of P limitation is caused by the low P sorption capacity of temperate and boreal soils.This can be seen in the simulation where all soils have high sorption capacities (CP ox ).In this simulation P limitation is comparable in strength to N limitation (Fig. 5).The shift is also robust against the assumption of increasing N deposition rates in the future.When N deposition rates are maintained at the present day level during the whole simulation (CN dp ), N limitation is still less than P limitation in the CP simulation (Fig. 5).We do not account for the increase in substrate availability for mineralization due to the thawing of boreal soils (Anisimov et al., 1997).This may partly counteract a shortage of nutrients at latitudes above 50 • N.
In summary, our findings suggest that high latitude ecosystems could become significantly more limited by both N and P in the future as growing seasons expand and soil stocks increase.From a theoretical point of view, a shift from N to P limitation is plausible as the buffer of available P, P adsorbed, gets depleted.However, it is still uncertain how soil C stocks will be affected by climate change, as the temperature dependence of heterotrophic respiration, which is also the main control of nutrient mineralization in JSBACH, is still a major uncertainty in simulating the C and nutrient cycles (Raddatz et al., 2007;Reich, 2010).

Conclusions
The quantification of P limitation remains challenging because of limited data and process understanding.Our simplistic approach does not account for all synergies between the element cycles.Using data from 34 separate fertilization studies, Marklein and Houlton (2012) found that an increase in N availability may accelerate P cycling via a stimulation of biochemical mineralization.Since the effect has not yet been quantified and the responses varied between plant and microbial mediated mineralization rates, the incorporation of such an interaction in a global model is not feasible at the moment.Moreover, it is commonly assumed that N fixation may be controlled by P availability (Cleveland et al., 1999;Wang et al., 2010).We use an empirical N fixation model which does not represent effects of changing P availability on N fixation, because the regulation mechanisms of N fixation are elusive and recent findings from a long term fertilization experiment reveal that N fixation in an intact tropical forest responded to molybdenum but not to P addition (Barron et al., 2009).Recently, Morford et al. (2010) raised the possibility that bedrock N input, which is omitted in the model, may present an important component of the N cycle.However, it is yet unclear whether this finding can be generalized from two forest sites to global scale ecosystems.Due to our limitation concept, additional nutrients in general increase productivity just in the case of CO 2 -induced nutrient limitation, which is relatively minor at present conditions.Overall, our model may tend to maximize nutrient limitation, and actual nutrient limitation may be less than projected.
We identified two processes which critically determine the strength and, to a lesser extent, the global distribution of P limitation: biochemical mineralization and the sorption capacity of soils.Both processes control the shift of P from unavailable forms (P adsorbed and in soil organic matter) to available ones (labile P), and therefore counteract a scarcity of P resulting from an accumulation of P in vegetation and soil organic matter.As both processes are poorly constrained by measurements (Coss and Schlesinger, 1995;Johnson et al., 2003) and changes in soil properties like maximum sorption capacity can not be captured by the empirical sorption model, more basic research on these soil processes is needed to better constrain P cycling.
Even so, our results suggest that models that do not account for P limitation overestimate the land C uptake during the 21st century.The geographic pattern of N and P limitation is to a large extent robust against the assumptions about soil P sorption, biochemical mineralization, CNP ratios, and N deposition.The reduction of tropical NPP by P limitation is in all cases larger than the reduction by N, while a shift from N to P limitation may occur with warming at higher latitudes.The low soil P sorption capacity of temperate soils is responsible for this shift in limitations.
A reduction in the response of NPP to increasing atmospheric CO 2 concentration due to P limitation is expected to weaken the land carbon sink.As a consequence the climate-C cycle feedback (Friedlingstein et al., 2006) increases, leading to a more pronounced warming than projected by C-only models.As most climate C cycle models do not account for P limitation, they most probably overestimate tropical C storage potential.This is problematic because tropical regions dominate the exchange of C between land and atmosphere in observation (Pan et al., 2011) as well as in model projections for the 21st century (Raddatz et al., 2007;Hickler et al., 2008).Therefore, it is important to include the P cycle in global models used for C cycle projections.of parameter samples, thus avoiding impracticable computational efforts due to numerous model simulations.We sampled the C : N ratios of the active, wood and slow pool, as well as the N : P ratios of the active and the wood pool.The N : P ratio of the slow pool was excluded from sampling, because this parameter affects the calibration of the biochemical mineralization and the calibration procedure is computationally too demanding to be performed.The stoichiometric parameters of the litter pools were not sampled, but the fraction of nutrients retranslocated prior to loss were.This was done to avoid correlations between parameters, as the C : nutrient ratio of litter must be larger than the ratio of the corresponding vegetation pool.For all 7 parameters, the probability functions were assumed to be uniform, and the ranges were constrained from literature (Table C1).Due to the scarce data, we assumed that the ranges of N : P ratios for wood and woody litter are comparable to the range of N : P ratios for leaves and litter, respectively.The retranslocation fractions were chosen such that the resulting litter   (2004) stoichiometries are in the range of observations (Table C1) (0.45-0.55 for the active pool and 0.5 for the wood pool).A set of 140 simulations was performed for each of the three nutrient setups (more simulations were not possible due to the high computational effort).All simulations were performed in a similar manner as the standard simulations.
The experiments showed that the overlaps between the CN, CP and CNP simulations are marginal (Fig. C2).The ranges of simulations within the 10th and the 90th percentiles do not overlap at all.In summary, the results support our conclusion that the reduction of land C uptake due to nutrient limitation is robust regarding the stoichiometric parameterization.

Fig. 1 .
Fig. 1.Schematic representation of pools and fluxes in JSBACH.Solid arrows indicate carbon fluxes and dashed arrows nutrient fluxes.The plant compartment consists of the three C pools: active (leaves and non-lignified tissue), wood (stem and branches) and reserve (sugar and starches).The litter compartment consists of non-lignified litter, and woody litter (lignified litter and fast-decomposing soil organic matter).The soil compartment consists of one pool (slow) representing slow-decomposing organic matter.All carbon pools except the reserve pool have a corresponding nutrient pool with variable C : N : P ratio (slow, non-lignified litter) or fixed C : N : P ratio (rest).There is one mobile plant pool representing mobile nutrients stored internally in plants.Soil mineral nitrogen is represented by a single pool (soil mineral pool), while mineral P is represented by labile (available) pool and sorbed pool.The opposing triangles indicated that the flux is controlled by phosphorus (red triangles), nitrogen (blue triangles) or both (mixed triangles).

Fig. 2 .
Fig. 2. The simulated P in vegetation (top left), sorbed P (top right), P in soil organic matter and litter (bottom left), and the annual P uptake by vegetation (bottom right) for present day.

Fig. 3 .
Fig. 3.The simulated change in land carbon storage under the SRES A1B scenario.Shown are the 10-yr mean of soil temperature (a), the CO 2 concentration as used in the forcing simulation (b), the resulting change in total land C storage (c), and the changes in the two main land compartments vegetation (d) and soil (e).

Fig. 4 .
Fig. 4. The reduction in C storage (kg m −2 ) by nutrient limitation at the end of the 21st century.Shown is the difference in the mean C storage (2070-2099) between the CN simulation and the C-only simulation (upper panel), and between the CP simulation and the C-only simulation (lower panel).The latitudinal means over land points are shown on the right side.

Fig. 5 .
Fig. 5. Mean reduction of NPP (%) due to CNL (N or P) averaged over 30 yr (2070-2099 (left) and 2320-2349 (right)) as latitudinal means.Solid line are results from the CN and the CP simulation for N and P limitation, respectively.The dashed lines are results from the sensitivity experiments CP ox , CP mo , and CN dp .

Fig. C1 .
Fig. C1.The simulated change in land carbon storage under the SRES A1B scenario.Shown are the 10-yr mean of soil temperature (a), the CO 2 concentration as used in the forcing simulation (b), the resulting change in total land C storage (c), and the changes in the two main land compartments, vegetation (d) and soil (e).CN-FLX and CP-FLX are the new simulations with flexible plant stoichiometry.

Fig. C2 .
Fig. C2.The simulated change in land carbon storage in the parameter pertubation simulations.Shown are the median (bold black), the range between the 10th and 90th percentile (color), and the range between mininum and maximum for each ensemble (grey).

Table C1 .
Ranges used to constrain the probability functions of the stoichiometric parameters in the LHS.
W P released from primary minerals by weathering [mol (P) m −2 s −1 ] depends on soil order Wang et al. (2010) M min minimal (background) biochemical mineralization rate [mol (P) m −2 s −1 ] depends on soil order (calibrated)

Table 3 .
List of experiments performed.The climatic forcing is the same for all simulations.See Appendix C for further simulations regarding the stoichiometric parameterization.

Table 5 .
Simulated P stocks (Gt (P)) from the simulation CNP compared with published estimates.The simulated stock sizes are means of the years 1970-1999.The litter pool contains leaf litter and fast decomposing soil organic matter.