A representation of the phosphorus cycle for ORCHIDEE (revision 4520)

. Land surface models rarely incorporate the terrestrial

The model captures the observed differences in the foliage stoichiometry of vegetation between an early (300-year) and a late (4.1 Myr) stage of soil development.The contrasting sensitivities of net primary productivity to the addition of either nitrogen, phosphorus, or both among sites are in general reproduced by the model.As observed, the model simulates a preferential stimulation of leaf level productivity when nitrogen stress is alleviated, while leaf level productivity and leaf area index are stimulated equally when phosphorus stress is alleviated.The nutrient use efficiencies in the model are lower than observed primarily due to biases in the nutrient content and turnover of woody biomass.
We conclude that ORCHIDEE is able to reproduce the shift from nitrogen to phosphorus limited net primary productivity along the soil development chronosequence, as well as the contrasting responses of net primary productivity to nutrient addition.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Introduction
As it has been acknowledged that human activity is changing Earth's climate, it is argued that climate research needs to sharpen its view (Marotzke et al., 2017).As its new suggested focus is the fate of the emitted carbon, which is closely linked to the extensive scientific debate about the importance of nutrient limitation (nitrogen and phosphorus supply) for future land carbon uptake (for example, Peñuelas et al., 2013;Wieder et al., 2015;Brovkin and Goll, 2015).Yet none of the Earth system models (ESMs), which are major tools in advancing the understanding of the role of human activities in the climate system, incorporate a terrestrial phosphorus cycle.
The few existing land surface models (LSMs) which account for interactions between phosphorus availability and the land carbon cycle suggest a significant role of phosphorus availability for ecosystems on highly weathered soils (Wang et al., 2010;Goll et al., 2012;Yang et al., 2014) with increasing significance as carbon dioxide concentration rises (Goll et al., 2012).However, these findings are highly uncertain due to processes which are poorly constrained by current observational data: soil phosphorus sorption dynamics, phosphatase-mediated mineralization, stoichiometric plasticity, leaf nutrient recycling, and the effects of phosphorus limitation on vegetation (photosynthesis, growth, allocation, mortality) (Wang et al., 2010;Goll et al., 2012;Yang et al., 2014;Reed et al., 2015).
Ecosystem manipulation experiments are shown to provide useful information to assess and evaluate LSMs (Medlyn et al., 2015;Meyerholt and Zaehle, 2015), which, in return, facilitate the interpretation of observation data and can guide the design of experiments (Medlyn et al., 2015).The long-term (6-10-year) fertilization experiment in a soil formation chronosequence in Hawaii (Harrington et al., 2001;Ostertag, 2001), with its contrasting availabilities of nitrogen and phosphorus along a soil age gradient going from young phosphorus-rich and nitrogen-poor soils to old highly weathered soils low in phosphorus but rich in nitrogen, provides an ideal test case for the evaluation of nutrient components in LSMs (Wang et al., 2007;Yang et al., 2014).
The potentially important influence of phosphorus availability on the land carbon balance and the recently initiated ecosystem-scale manipulation experiments in phosphoruspoor environments (for example, Ama, 2017;Euc, 2017;IMB, 2017;AFE, 2017) as well as other projects related to the role of phosphorus in ecosystem functioning (for example, SPP, 2017;QUI, 2017), call for the need for new phosphorus enabled LSMs to keep track of these actions (Reed et al., 2015).
Here, we describe the implementation of the terrestrial phosphorus cycle in the ORCHIDEE LSM (Krinner et al., 2005) following the principles developed earlier for the introduction of the nitrogen cycle into ORCHIDEE (Zaehle and Friend, 2010).It is the first global phosphorus model which explicitly simulates root uptake of dissolved phosphorus accounting for soil moisture effects on soil phosphorus mobility.The model is then evaluated with data from a long-term fertilization experiment in a soil formation chronosequence in Hawaii (Harrington et al., 2001;Ostertag, 2001).

Model description
The land surface model used for this study, ORCHIDEE, is based on two different modules (Krinner et al., 2005, their Fig. 2).The first module describes the fast processes such as the soil water budget and the exchanges of energy, water, and CO 2 through photosynthesis between the atmosphere and the biosphere.The second module simulates the carbon dynamics of the terrestrial biosphere and essentially represents processes such as maintenance and growth respiration, carbon allocation, litter decomposition, soil carbon dynamics, and phenology.Global vegetation is described by 13 metaclasses which correspond to plant functional types (PFTs) with a specific parameter set (one for bare soil, eight for forests, two for grasslands, and two for croplands).
The major modifications since Krinner et al. (2005) are listed in the following: a slightly revised carbon allocation scheme from a recent side branch of ORCHIDEE (Naudts et al., 2015), which avoids the capping of the leaf area index at a predefined value; an explicit representation of mesophylic conductance to CO 2 and omission of direct effects of soil moisture stress on the maximum rate of carboxylation (V cmax ) (Vuichard, unpublished data); and a revised thermodynamic scheme which accounts for the heat transported by liquid water into the soil, in addition to the heat conduction process (Wang et al., 2016).

Starting version
The implementation of the phosphorus (P) cycle in OR-CHIDEE was done in the nitrogen enabled version of OR-CHIDEE (ORCHIDEE-CN) (Vuichard, 2017 unpublished data).ORCHIDEE-CN is a re-implementation of the nitrogen cycle from a discontinued version of ORCHIDEE (which became OCN, Zaehle and Friend, 2010;Zaehle et al., 2011) in a recent version of ORCHIDEE (r3623).The nitrogen cycle in OCN is well evaluated (De Kauwe et al., 2014;Zaehle et al., 2014;Walker et al., 2015;Meyerholt and Zaehle, 2015) and identical to the one in ORCHIDEE-CN except for the parametrization of the relationship between leaf nitrogen concentration and maximum carboxylation capacity of photosynthesis (V cmax ) as ORCHIDEE (r3623) uses a different carbon assimilation scheme than originally used in Zaehle and Friend (2010).V cmax is directly derived from the leaf nitrogen concentration at the respective canopy level following Kattge et al. (2009): Adsorbed labile phosphorus can be transformed into more recalcitrant forms which are irreversibly lost to biota.When plants take up phosphorus from soils it enters the plant labile phosphorus pool, from which it is allocated to growing plant tissues.The plant labile phosphorus concentration is buffered by a reserve pool which serves as a long-term storage to buffer seasonal variations in phosphorus demand and supply.When plant tissue is shed, part of the phosphorus is recycled (broken lines), while the rest enters the litter pools, from where it is either transformed into soil organic matter or mineralized.
where N * leaf,h is nitrogen concentration in leaves at canopy level (h).N * leaf,h is derived from nitrogen in leaf biomass per ground area (N leaf ) using an exponential canopy nitrogen profile (Johnson and Thornley, 1984).N * leaf,h is corrected for a certain fraction of structural nitrogen per leaf carbon (N str ) which does not contribute directly to the carboxylation capacity of photosynthesis (Table 2): where κ N is the extinction coefficient.The electron transport capacity (J max, h ) is derived from V cmax,h using the relationship from Kattge and Knorr (2007) which accounts for acclimation of photosynthesis to monthly temperatures.
In the following the representation of the terrestrial phosphorus cycle and its interaction with the cycles of carbon and nitrogen are described.All variables and parameters can be found in Tables 1 and 2, respectively.

Vegetation: phosphorus uptake, allocation, and turnover
Vegetation biomass (carbon, nitrogen, and phosphorus) is separated into leaves, roots, sapwood, heartwood, and short-term (labile) and long-term storage (reserves) (Zaehle and Friend, 2010;Naudts et al., 2015) (Fig. 1).We prescribed boundaries for the stoichiometry of leaves, roots, sapwood, and heartwood, but not for labile and reserves (Table A1) (Zaehle and Friend, 2010).Plants take up labile phosphorus dissolved in soil solution (solution phosphorus) following the representation of root nitrogen uptake (Zaehle and Friend, 2010).Plant nutrient uptake is simulated as a function of mineral nutrient availability with the aim to account for the increase in uptake in nutrient starved plants by increasing the uptake capacity per root surface (Schachtman et al., 1998), as well as indirectly through increased root growth and exploration of the soil by roots to increase their resource acquisition (Schachtman et al., 1998).Mycorrhizae are implicitly included in root biomass as mycorrhizal hyphae show comparable uptake characteristics as roots (Schachtman et al., 1998).As the concentration of phosphorus in roots is orders of magnitudes larger than the concentration in the soil solution, passive uptake of phosphorus via diffusion is negligible (Schachtman et al., 1998).Thus, only active uptake via specialized transporters on the root surface is accounted for in the model.Hereby, the model does not distinguish between organic and inorganic forms of dissolved labile phosphorus.
www.geosci-model-dev.net/10/3745/2017/Geosci.Model Dev., 10, 3745-3770, 2017 The x root uptake capacity (u max ) per root mass (C root ) for a given solute concentration follows the combined behavior of low-affinity and high-affinity transporter systems working in parallel, which typically shows no saturation at high soil solute concentrations (Kronzucker et al., 1995;Zhang et al., 2009) and is given by where v max is the maximum uptake capacity, a root P sol is the dissolved labile phosphorus concentration at the root surface (Eq.23), and c k a unit conversion factor using the soil-typespecific parameter for soil moisture content at saturation in ORCHIDEE (mcs) as an approximation of pore space following Smith et al. (2014).
The combined behavior of the two uptake systems is approximated by the term in brackets, where the linear factor (k P min ) was chosen to match the observed rate of increase in overall phosphorus uptake at high dissolved labile phosphorus concentration (low-affinity transporter) (Zhang et al., 2009) (Table 2).The values for the Michaelis-Menten constants are averages of the values reported in Schachtman et al.
(1998) (page 448) for k P min of the high-affinity system (Table 2).We initially used the values reported in Bouma et al. (2001) for v max for orange trees, but had to reduce these values by a factor of 10 to achieve realistic uptake behavior (Table 2).Plant uptake (F up ) is derived from multiplying the root uptake capacity by the root carbon mass (C root ) and is scaled with f PN plant to account for actual phosphorus demand and with f temp to avoid phosphorus accumulation in plants and soil at low temperature: The temperature scaling function (Zaehle and Friend, 2010) is given by As phosphorus uptake is energetically costly (Schachtman et al., 1998), plant phosphorus uptake is down-regulated ac-cording to the P : N ratio of plant tissue (pn plant ), avoiding excessive uptake of phosphorus (luxury consumption) when tissue phosphorus concentrations are at the prescribed maximum (pn leaf, max ): The dependency of phosphorus uptake on pn plant is described as pn leaf, min − pn leaf, max , 0.0 , 1.0 (7) where pn leaf, min and pn leaf, max are the minimum and maximum foliage P : N ratios.Maximum uptake is reached when pn plant equals pn leaf, min by the use of a minimum function.
Note that because neither the nitrogen nor the phosphorus concentration in the plant labile phosphorus pool (P labile ) is constrained by a prescribed P : N ratio, the actual value of f PN plant may be higher than 1.Further, we scale plant phosphorus uptake with a temperature function (f temp ).We use the same equation as is used to www.geosci-model-dev.net/10/3745/2017/Geosci.Model Dev., 10, 3745-3770, 2017 Figure 2. The effects of nitrogen and phosphorus stress (sub-optimal internal availability) on vegetation and associated processes in OR-CHIDEE.A shortage of internal nutrients reduces tissue nutrient concentrations, overall growth (Eq.12), the shoot-to-root ratio of new growth (Eq.14), litter quality, and, in the case of nitrogen, the carboxylation efficiency of photosynthesis (Eq.1).In addition, processes enhancing the availabilities of nitrogen and phosphorus are up-regulated (Eqs.18 and B2).
scale soil carbon decomposition (Krinner et al., 2005), phosphorus mineralization, biochemical mineralization (Eq.18), and nitrogen uptake and mineralization (Zaehle and Friend, 2010).The equation avoids the accumulation of phosphorus in plants or soils at low temperatures, mimicking the inhibition of biological processes when soils are frozen, which is not explicitly represented in ORCHIDEE.
The phosphorus taken up by plants enters their labile phosphorus pool (P labile ) whose dynamics are given by where τ i is the fraction of foliage or roots shed each time step, f trans,i is the fraction of phosphorus recycled and transferred to plant labile phosphorus before tissue is shed, G P is labile phosphorus allocated to new biomass, and F reserve is the flux to or from the long-term storage (P reserve ).
Following the dynamics of labile nitrogen (Zaehle and Friend, 2010), P labile is limited to a maximum size P labile, max which is taken as the phosphorus required to allocate the entire labile carbon pool according to the current growth rate and the maximum foliage phosphorus concentration.Any excess labile phosphorus is transferred to P reserve and is mobilized again if the size of the labile phosphorus pool falls below P labile, max : Following the assumption regarding nitrogen concentration (Zaehle and Friend, 2010), the phosphorus concentration in newly grown plant tissue is assumed to depend directly on the phosphorus concentration in the plant labile pool, providing a link between tissue activity and plant labile phosphorus availability.Foliage phosphorus concentration changes are simulated explicitly, with the phosphorus content of nonfoliage tissue varying in proportion to that of the foliage, as observed along gradients of soil fertility (Heineman et al., 2016).The phosphorus required (G P ) to sustain the current growth (G C ) of new tissue can therefore be written as where f i are the fractions of carbon allocated to foliage (i = leaf), roots (i = root), fruits (i = fruit), and sapwood (stalks for grass) (i = sap) which are calculated dynamically (Zaehle and Friend, 2010;Naudts et al., 2015), nc leaf and pn leaf are the nitrogen to carbon and the phosphorus to nitrogen ratio of current foliage, λ i are the phosphorus to carbon allocation to tissue i relative to the phosphorus to carbon allocation to leaves (λ leaf = 1), and D leaf is an empirical elasticity parameter.Analogous to leaf nitrogen concentrations (Zaehle and Friend, 2010), the foliage phosphorus concentrations are dynamic state variables.If the plant labile phosphorus pool is not sufficient to maintain the current phosphorus concentration at the current carbon growth rate G C , the phosphorus concentration of newly grown leaf tissue is allowed to decrease relative to the concentration of existing foliage.Conversely, if plant labile phosphorus is larger than required, and the plant is not in the phase of flushing new foliage, P concentrations are allowed to increase.
To dampen day-to-day variations in tissue nutrient concentrations, such as at the beginning of the growing season, an empirical elasticity parameter (D leaf ) is included.Meyerholt and Zaehle (2015) tested different assumptions about the stoichiometric flexibility in the OCN model and showed Geosci.Model Dev., 10, 3745-3770, 2017 www.geosci-model-dev.net/10/3745/2017/ that stoichiometric flexibility is needed to reproduce observational data from fertilization experiments; however, they found that the original formulation in OCN was too flexible.Vuichard et al. (unpublished data) revised the formulation of the dampening equation for the leaf nitrogen concentration which is also applied for the leaf phosphorus concentration: We adapted the dependency of biomass growth on plant labile nitrogen availability (Zaehle and Friend, 2010) for the dependency on plant labile phosphorus availability (Fig. 2): if the plant labile phosphorus concentration in vegetation fails to match the P requirement of biomass carbon growth, the growth of plant tissue is reduced proportionally to match phosphorus availability in the plant labile pool: where g max is a unit-less scalar regulating the maximal daily fraction of P labile allocated to growth, to avoid a complete depletion of P labile at any given time step.g max is also used to regulate the allocatable fraction of labile carbon and labile nitrogen, and it is a function of temperature (Naudts et al., 2015).G P is the estimated amount of phosphorus needed to support growth.C growth (G C ) is then scaled by the minimum of growth limitation factors derived from phosphorus availability (p lim ) and nitrogen availability (n lim ) (see Eq. 22 in the Supplement of Zaehle and Friend, 2010): Nutrient stress in general affects the ratio leaf to root portioning of new growth: where f LF is a function relating leaf mass to root mass based on the pipe theory (Shinozaki et al., 1964) as originally implemented by Zaehle and Friend (2010) and recently updated by Naudts et al. (2015).n scal is the actual nutrient stress factor and is derived from the minimum of the nitrogen (n scal,N ) and phosphorus (n scal,P ) stress scaling factors: n scal,P is given by the deviation of the actual plant phosphorus concentration from the maximal leaf phosphorus concentration relative to carbon concentration: where pc leaf, ave is the average of the maximum and minimum leaf phosphorus to carbon ratios (pn leaf, min , pn leaf, max , nc leaf, min , and nc leaf, max ) and pc plant the growing season average of the plant labile phosphorus to labile carbon concentration: The calculation of n scal,N follows the calculation of n scal,P with the exception that the deviation of the actual plant nitrogen concentration from the maximal leaf nitrogen concentration relative to carbon concentration is used (Zaehle and Friend, 2010).Turnover of biomass phosphorus follows strictly the turnover of each biomass pool as described in Krinner et al. (2005).The phosphorus fluxes are derived from the carbon fluxes and the corresponding stoichiometric ratios, subtracting a fixed fraction of the phosphorus which is resorbed and added to the plant labile pool.

Litter and soil organic matter
The turnover of litter and soil organic matter follows the CENTURY model (Parton et al., 1993), which describes decomposition as a function of substrate availability, clay content, soil moisture, and soil temperature.Organic matter is separated into structural and metabolic litter and three soil organic matter pools (fast, slow, passive) which differ in their respective turnover times with no vertical discretization.Due to the fast turnover of microbial communities, microbial biomass is assumed to be always adjusted to the availability of labile organic matter and is thus part of the fast soil organic matter pool.The model is described in detail elsewhere (Krinner et al., 2005;Zaehle and Friend, 2010).The nitrogen concentrations of decomposing material are assumed to vary linearly with soil mineral nitrogen content.Instead of applying a comparable (empirical) approach for the phosphorus concentration of decomposing material (Parton et al., 1993;Kirschbaum et al., 2003), the phosphorus concentrations vary mechanistically as a function of biochemical mineralization (Eq.18) (Wang et al., 2010;Goll et al., 2012;Yang et al., 2014).

Biochemical mineralization
Biochemical mineralization (phosphatase-mediated) decouples the mineralization of phosphorus partly from carbon decomposition and nitrogen mineralization (McGill and Cole, 1981).In contrast to biological mineralization of nitrogen and phosphorus, biochemical mineralization is not driven by the energy demand of microorganisms.Although phosphatase activity, which is a qualitative measure of biochemical mineralization, is common in soils (Stewart and Tiessen, 1987), the quantification of the mineralization rates in the field is not yet possible.
www.geosci-model-dev.net/10/3745/2017/Geosci.Model Dev., 10, 3745-3770, 2017 We simulate biochemical (phosphatase-mediated) mineralization of phosphorus (F bcm ) with the aim to account for the observed increase in F bcm when plants experience suboptimal P-to-N availabilities as an approximation of the stoichiometric status of the whole ecosystem (Margalef et al., 2017), including the effect of substrate availability on mineralization (McGill and Cole, 1981).
where τ x,ref is the turnover time of phosphorus in soil organic matter pool x (P x ) and the other two variables are scaling functions.First, biochemical mineralization is scaled according to the P : N status of vegetation (f PN plant ) to account for the observed link between rizosphere phosphatase activity and plant nutritional status (Fox, 1992;Hofmann et al., 2016).Second, it is scaled with the same equation used to scale mineralization and root uptake according to soil temperature.The values of τ x,ref are set arbitrarily, due to the lack of observational constraints, to half the turnover times used for the "biological mineralization" of organic matter (Krinner et al., 2005), except for τ passive, ref , which is set to zero to account for inaccessible phosphorus in stabilized nutrientrich organic matter (Tipping et al., 2016).

Soil mineral phosphorus
The release of phosphorus from primary minerals is the primary source of phosphorus for many terrestrial ecosystems.In this study, we prescribed site-specific release rates (F weath ), but a dynamic phosphorus weathering routine is implemented in ORCHIDEE which is described in the Appendix.
Leaching (F leach ) of dissolved phosphorus (P sol ) occurs in proportion to the fraction of soil water ( ) lost by the sum of simulated drainage and surface runoff (q): We assume that at each time step a fixed fraction (k s ) of dissolved labile phosphorus is adsorbed onto soil particles and the remaining fraction (1 − k s ) is dissolved.Instead of the commonly used Langmuir equation, we chose a linear approach for sorption, which works well for low soil phosphorus concentrations, which are common in most natural ecosystems (McGechan and Lewis, 2002).The calibration of the Langmuir equation for global application represents a major challenge as global datasets on soil phosphorus content are limited (Yang and Post, 2011) and parameters cannot be derived with enough confidence.Given the high sensitivity of the dynamics of available phosphorus on the sorption dynamics, we choose a simple but sufficiently constrainable approach.Thus the dynamics of sorbed labile phosphorus (P sorb ) are given by The dynamics of dissolved labile phosphorus (P sol ) are given by (see Appendix A in Goll et al., 2012, for details) where F weath (Eq.A1) is phosphorus release from primary minerals, F up (Eq.4) is phosphorus uptake by plants, F bcm (Eq.18) is biochemical mineralization, and F min is the biological mineralization of phosphorus (Parton et al., 1993).The fraction of adsorbed to total soil labile phosphorus is derived from a global dataset of soil phosphorus fractions (Yang et al., 2013) and we use USDA soil-order-specific parameter values.Further, we assume a constant rate at which adsorbed mineral phosphorus becomes strongly sorbed (τ sorb ) and is subsequently fixed into secondary minerals.The turnover time of sorbed phosphorus with respect to occlusion is derived from the difference in occluded phosphorus among the sites of the Hawaii chronosequence (Violette, unublished data).

Root zone mineral phosphorus
As the mobility of phosphorus in soil is very low, plant uptake tends to be limited by the replenishment of phosphorus to the root surface rather than by the root uptake capacity itself (Schachtman et al., 1998).We simulate the labile phosphorus concentration in soil solution in root contact as a function of plant uptake and diffusion of phosphorus from the surroundings towards the root surface without a vertical discretization.We assume that plant uptake is small compared to the actual amount of dissolved phosphorus in total soil volume (Johnson et al., 2003), and thus its effect on the dissolved phosphorus concentration is limited to a small band around the surface of roots.The diffusion of phosphorus from the surroundings to the root surface (F diff ) follows Fick's law: where D is the permeability of the soil to phosphorus and P sol is the difference in the phosphorus concentrations between the soil solution at the root surface (a root P sol ) and the solution in the surrounding soil volume outside the diffusive zone around the root ( P sol ).Assuming a homogeneous distribution of soil water, changes in the phosphorus concentration in the root zone are given by P sol = (a root − 1) P sol (23) where is the volumetric soil water content and a root is the relative reduction of labile phosphorus in soil solution at the root surface compared to the surroundings.As a root ≤ 1, the diffusion is a single direction flux.
Geosci.Model Dev., 10, 3745-3770, 2017 www.geosci-model-dev.net/10/3745/2017/ The permeability D is calculated analogously to the diffusion coefficient of phosphorus in soils following Barraclough and Tinker (1981), which accounts for the increased path length in soil using a tortuosity factor (tf). D is further corrected for half the path length between uniform distributed root cylinders (r diff ): where D 0 is the diffusion coefficient in free water, the volumetric soil water content, and c a unit conversion factor.The tortuosity factor is given by a broken-line function of (Barraclough and Tinker, 1981) where l is the soil water content at which the two functions intersect, and f 1 and f 2 are empirical parameters (Barraclough and Tinker, 1981).
We assume that the diffusion path (r diff ) can be approximated by half the distance between uniformly distributed roots.We restrict the diffusion path length to 0.1 m, as the effect of active root phosphorus uptake on the soil phosphorus concentration at a distance of more than 10 cm is negligible (Li et al., 1991).Following Bonan et al. (2014), we derive half the distance between roots as r diff = min 0.1, (π RLD) 0.5  (26)   where the root length density (RLD) (root length per volume of soil) is given by where r d is the root-specific density and π r 2 r is the crosssectional area calculated from the fine root radius, r r , and M * root is the root biomass density in the soil volume.
The change in the difference in the dissolved labile phosphorus concentration between the root surface and the surroundings, a root , is then derived by where F up is plant uptake of phosphorus as described earlier.

Competition between microbes and plants
The competition between microbes and plants for dissolved labile phosphorus is handled analogously to the competition for soil mineral nitrogen (Zaehle and Friend, 2010).Gross phosphorus immobilization, gross biological phosphorus mineralization, biochemical mineralization, as well as plant phosphorus uptake are calculated half-hourly.At any time step, immobilization due to litter and soil organic matter decomposition is given priority in accessing nutrients from gross biological mineralization.This is in line with recent findings regarding the variability in the nitrogen use efficiency of microbes (Mooshammer et al., 2014), which indicates a dominance of microbes in accessing soil nitrogen and results in increasing immobilization with decreasing litter nutrient content.
The nutrient requirement for the build-up of soil organic matter, which affects the nutrients retained from litter decomposition, is dependent on the C : N : P ratio of soil organic matter, whereas the C : N ratios depend on the soil mineral nitrogen concentration (Zaehle and Friend, 2010).Increasing plant uptake of nitrogen reduces the soil mineral nitrogen concentration and thereby reduces the nitrogen retained from litter decomposition in soil organic matter due to its effect on soil C : N ratios.

Input fields
The parameter describing soil phosphorus sorption (k s ) is USDA soil order specific.The parameters for phosphorus release from minerals (s shield , w l , E a, l ) are lithological class specific and read in from the GliM lithological map (Hartmann and Moosdorf, 2012).

Site-scale simulation
The long-term field fertilization experiment along the Hawaiian soil development chronosequence provides an ideal test case for the nutrient components of ORCHIDEE (Vitousek, 2004).We selected sites for which sufficient observational data are available (Harrington et al., 2001;Ostertag, 2001): a 300-year old site which is nitrogen limited (Thurston) and a 4.1 Myr old site which is phosphorus limited (Kokee).The two sites have similar climatic conditions (Table 3) and are dominated by the same tree: Metrosideros polymorpha (Crews et al., 1995).
We run the model with observed meteorological data (Harris et al., 2014;Thornton et al., 2016) prescribing nutrient boundary conditions, namely inputs of phosphorus and nitrogen by atmospheric deposition (Chadwick et al., 1999) and inputs of phosphorus by weathering (Crews et al., 1995).To do so, we deactivated the module for dynamic phosphorus weathering (see the Appendix) and instead prescribed a constant site-specific release rate.In addition, we prescribe sitespecific physico-chemical soil properties (Crews et al., 1995;Chorover et al., 2004;Olander and Vitousek, 2004).The prescribed vegetation cover is tropical evergreen broadleaf vegetation.
We equilibrated the biogeochemical cycles of the 4.1 Myr old site to the climatic conditions and the nutrient inputs using the semi-analytical spinup procedure (Naudts et al., 2015) which was extended to handle nutrient cycles.To capwww.geosci-model-dev.net/10/3745/2017/Geosci.Model Dev., 10, 3745-3770, 2017 5.9 5.9 Harrington et al. (2001) ture the transitional nature of the 300-year old site, we perform a 230-year long spinup simulation.The differences between the simulation duration and the actual age of the site are due to a correction for an initial amount of biomass we have to set in ORCHIDEE for technical reasons (see Appendix C).We extended the spinup simulations of both sites into a set of three nutrient addition simulations: adding only nitrogen, only phosphorus, and nitrogen and phosphorus together.A total of 10 g (N) m −2 yr −1 and 10 g (P) m −2 yr −1 are added in the model simulations homogeneously distributed across the year.In the field the same amount per year was applied, but semi-annually (Harrington et al., 2001;Ostertag, 2001).

Forcing data
The meteorological forcing data for ORCHIDEE are derived from Daily Surface Weather Data on a 1 km Grid for North America (DAYMET), version 3 (Thornton et al., 2016).The data include meteorological information (short-wave radiation, maximum daily temperature, minimum daily temperature, daily precipitation sum) for the period 1980-2013.In DAYMET, the mean annual surface temperature and annual sum of precipitation for the actual locations of the two sites (Table 3) substantially deviate from the values reported at the sites (Crews et al., 1995).Therefore, we pick the closest site nearby, Thurston (16 km distance; lat = 19.8318,lon = −155.411),which has an annual sum of precipitation of 2500 ± 250 mm yr −1 and an average annual temperature 16 ± 1 • C as reported in Crews et al. (1995), and use it for both sites (as the DAYMET data for the Kokee island did not include any grid point with an appropriate climate).
We extract additional information which is needed to run ORCHIDEE, namely surface pressure, long-wave downward radiation, and wind, from a 0.5 × 0.5 • reanalysis dataset (CRU-NCEP) (Harris et al., 2014) using the coordinates of the Thurston site.We correct surface pressure from CRU-NCEP with the actual altitude of each site using a lapse rate and derive specific humidity from water pressure, air temperature, and surface pressure.
The annual inputs of nutrients by atmospheric deposition and weathering (Table 3) are kept constant and are evenly distributed throughout the year.

Site-specific parametrization
We use site information collected from the literature (Harrington et al., 2001;Ostertag, 2001) to parametrize the model (Table 3).We account for differences in the soil characteristics between sites, but use a common parametrization for all biological processes.Thereby, we are able to evaluate the differences in vegetation due to differences in soil characteristics and chemical weathering solely.
As soils are not vertically discretized in ORCHIDEE, we average observations when given for different soil horizons.The soil fractions for Thurston are assumed to be 50 % sand, 25 % silt, and 25 % clay due to a lack of site-specific information.
We account for changes in the sorption characteristics of volcanic soils as they develop.For the 300-year site, we use the average value of k s for Andisols from Yang and Post (2011).For the 4.1 Myr site k s was scaled with the relative difference in soil phosphorus sorption capacity between the two sites as computed dynamically in the P-enabled ver-Geosci.Model Dev., 10, 3745-3770, 2017 www.geosci-model-dev.net/10/3745/2017/sion (Violette, unpublished data) of the WITCH mechanistic weathering model (Goddéris et al., 2006).P release rates from primary and secondary minerals are inferred from the observed differences in the chemical composition of minerals between sites along the chronosequence from Crews et al. (1995) (Violette, unpublished data).The reference decomposition rates of soil organic matter pools by biochemical mineralization (τ x,ref ) are chosen so that the nitrogen to phosphorus ratio of soil organic matter is close to the observation of approximately 10 g (N) g −1 ((P) for sites older than 10 kyr in the Hawaiian chronosequence (Crews et al., 1995).This is a common procedure (Wang et al., 2010;Goll et al., 2012;Yang et al., 2014) as this flux remains to be quantified in the field.
We prescribe observed values for specific leaf area, which is a fixed parameter and does not vary over time, and use the 25th and 75th percentiles of observed values of leaf P : N of the dominant tree species (Kattge et al., 2011) as boundaries for the leaf P : N ratio.We further increase the critical leaf age, which scales leaf turnover related to leaf age in ORCHIDEE, from 1.4 to 6 years to account for the substantially longer lifespan of leaves at both sites (Harrington et al., 2001) compared to the default value of ORCHIDEE (Naudts et al., 2015).
Following Yang et al. (2014), we adjust the turnover of the passive soil organic matter pools to achieve soil organic carbon stocks close to the observations.The same turnover rates are used for both sites.
The remaining parameters (including parameters for biological nitrogen fixation) are taken from the global parametrization of ORCHIDEE of the tropical evergreen broadleaf PFT (Table A1).

Analysis
We aggregate estimates of root production approximated by soil respiration from Ostertag (2001) and compare it to the simulated below-ground component of NPP (namely, NPP allocated to below-ground sapwood, below-ground heartwood, and fine roots).We sum the simulated above-ground component of NPP allocated to sapwood and heartwood and compare it to estimates of wood production based on wood increment and woody litter fall (Harrington et al., 2001).All other components of simulated NPP (fruit, leaf, reserve) are pooled and compared to estimates of non-woody NPP based on litterfall (Harrington et al., 2001).Simulated nutrient use efficiencies (NUE, PUE) are calculated as where NPP is annual NPP and X uptake the annual uptake of nutrient X, X ∈ {N, P}.Simulated nutrient use efficiencies are then compared to estimates derived from on-leaf litter fall, root growth, and wood increment in combination with the chemical composition of leaves and wood (Harring-ton et al., 2001).We further separated the nutrient use efficiencies into its underlying components, carbon production rate per biomass nutrient (X prod ) and nutrient residence time (X MRT ), following Finzi et al. (2007): where X content is the whole plant content (g (X) m −2 ) of nutrient X, X ∈ {N, P}.No information on below-ground productivity is available for the fertilization treatment (Harrington et al., 2001).Thus, we calculate above-ground nutrient use efficiencies (aNUE, aPNUE) using above-ground NPP instead of total NPP.Simulated apparent leaf lifespan (including climatic effects) is calculated by dividing the annual mean of leaf mass by the annual sum of NPP allocated to leaf growth.We calculate the rate by which soil organic phosphorus is biochemically mineralized and compare it to measurements of phosphatase activity (Olander and Vitousek, 2000) as a proxy of potential biochemical mineralization due to lack of alternatives.
The uncertainty ranges of simulated variables are given by the SD of annual values.We perform Student's t tests to determine whether the fertilization treatments resulted in significant differences in the tested variables compared to the control experiment in observations (when sufficient information is available) and simulations.

Control simulation
The comparison of carbon and nutrient cycle related ecosystem properties in the control simulations with observations allows us to detect model biases which facilitate the evaluation of the outcome of the fertilization experiments.
Net primary productivity (NPP) at both sites is well reproduced by the model (Fig. 3).At both sites NPP was not calibrated and thus is an independent model outcome.The simulated inter-annual variation in NPP at the 4.1 Myr site is more than twice as large as at the 300-year site.The model tends to capture the allocation pattern of NPP to the different plant tissues (Table 4).While wood growth is overestimated at both sites, the relative allocation to leaf and roots is rather well reproduced: the simulated ratios between root and leaf growth of 1.35 and 1.50 for the 300-year and 4.1 Myr sites, respectively, are close to the observed ratios of 1.36 and 1.33.This shows that the allocation scheme in ORCHIDEE, which accounts in a simplistic way for changes in the allocated fraction of NPP into below-ground allocation in response to stress (light, nutrient, water) (Zaehle and Friend, 2010;Naudts et al., 2015), gives reasonable results.
The simulated biomass stocks are in good agreement with the observations, with the exception of woody biomass at the 4.1 Myr site (Fig. 3).As wood growth is overestimated, the low woody biomass can be linked to an overestimation of wood turnover (Appendix D).Comparably, the slight overestimation of fine root biomass at the 300-year site is linked to an overly high turnover of fine roots (Appendix D).The large differences in observed fine root turnover between sites (Ostertag, 2001) cannot be captured by the model as fine root turnover is constant in ORCHIDEE.
Nutrient use efficiencies (NPP divided by plant nutrient uptake) are implicit plant properties that depend on the tissue stoichiometry, as well as the relative allocation of NPP to the various plant organs and their respective turnover rates.The nitrogen use efficiency (NUE) and phosphorus use efficiency (PUE) are underestimated between 20 and 50 % (Table 4).The analysis of the underlying components of nutrient use efficiencies following Finzi et al. (2007) (Eq. 30) indicates that the underestimation of NUE is mainly driven by the low carbon productivity per plant nitrogen (N prod ), while the low bias in PUE is due to a combination of low P prod and the short residence time of plant phosphorus (P MRT ) (Appendix D).The low N prod and P prod at the 300-year site can be attributed to the overestimation wood biomass and its nitrogen content (Appendix D).At the 4.1 Myr site the underestimation of the nutrient content of biomass and wood biomass has opposing effects on the nutrient use efficiencies.The general underestimation of the residence time of phosphorus (Table 4) is likely due to an underestimation of the phosphorus content of longlived plant tissue and the overestimation of wood turnover.Additionally, the extremely low concentration of plant available phosphorus at the 4.1 Myr site results in a set of physiological and morphological adaptation mechanisms which increase P MRT but are not resolved in ORCHIDEE (for example, changes in root morphology and turnover and leaf phosphorus recycling) (Schachtman et al., 1998;Niu et al., 2013;Reed et al., 2015).In ORCHIDEE the somewhat longer nutrient residences times at the young site are primarily due to the site's transient state in which biomass is still accumulating (Appendix D).
The simulated inter-annual variabilities in nutrient use efficiencies, in particular in the simulated NUE at the 300year site, are very large due to a substantial but highly variable contribution of nutrients from internal reserves to new biomass growth (not shown).Reserves can amount to up to 75 % of peak nutrient content in fine roots and leaves during Geosci.Model Dev., 10, 2017 www.geosci-model-dev.net/10/3745/2017/ the last growing season for evergreen plant functional types in ORCHIDEE (Zaehle and Friend, 2010).In the model, variations in reserves can be large in nutrient-poor environments which are subject to periods of reduced growth unrelated to nutrient starvation (here: drought).During droughts the nitrogen reserves are filled with foliage nitrogen which is recycled prior to leaf fall, while the nitrogen reserves are depleted when water availability and subsequent growth is high.Defoliation experiments indicate that plants can rely on substantial amounts of internally stored reserves of carbohydrate and nutrients which allow them to survive multiple defoliation events (Hartmann and Trumbore, 2016); however, the extent to which plants rely on internal storages is strongly species dependent (Piper and Fajardo, 2014), and the role of nutrients is often overlooked (Hartmann and Trumbore, 2016).Nonetheless, Ichie and Nakagawa (2013) showed that in Dryobalanops aromatica stored phosphorus accounted for 67.7 % of the total phosphorus requirements for reproduction, while stored N accounted for only 19.7 %, indicating substantial nutrient reserves in tropical trees.The simulated leaf N : P ratio is, at 15.2 at the 300-year site, somewhat higher than the observation of 12.6 ± 1.56, which lies just outside the lower boundary (np leaf, min ) of 12.83 in ORCHIDEE (Table 4).At the 4.1 Myr site, the simulated leaf N : P ratio is at 17.5 g (N) g −1 (P) very close to the observation of 17.3 g (N) g −1 (P).Foliage N : P ratios of less than 14 are commonly associated with nitrogen limitation and ratios above 16 with phosphorus limitation (Koerselman and Meuleman, 1996).Although the simulated leaf nitrogen concentration of 0.93 %(dryweight) at the 300-year site is 27 % higher than observed, it is substantially lower than the optimal concentration prescribed (cn leaf, min ) of 3.33 %dryweight, indicating substantial effects of the low soil nitrogen concentration (Table 4) on productivity and allocation at the 300-year site in the model.Thus, despite being useful tools, the common use of threshold stoichiometric ratios and models with rigid plant traits is somewhat limited when it comes to species-specific responses (Verheijen et al., 2016).
The differences in stoichiometry mirror differences in the respective availabilities of mineral nitrogen and soil labile phosphorus (Table 4).While the concentration of mineral nitrogen is extremely low at the young site due to a high immobilization demand of accumulating soil organic matter, the concentration is high at the old site, where immobilization demand is met by the mineralization of nitrogen from organic matter.In the case of phosphorus, the high phosphorus input of 434.0 mg m −2 yr −1 at the 300-year site keeps soil labile phosphorus concentration high despite the high immobilization demand.At the the old site, the extremely low phosphorus inputs of 0.27 mg m −2 yr −1 result in low soil labile phosphorus concentration as the ecosystem relies primarily on the mineralization of phosphorus from soil organic matter.
At the 300-year site, the simulated C : N : P stoichiometry of soil organic matter is at 309 : 16:1 nearly twice as rich in phosphorus as observed (425 : 28 : 1).As no significant differences in phosphatase activities among sites were observed (Ostertag, 2001) and the biochemical mineralization in OR-CHIDEE is calibrated to achieve realistic phosphorus concentration in soil organic matter in the long term (4.1 Myr site), the deviation of the simulated from observed phosphorus concentrations has processes other than biochemical mineralization.A recent data analysis suggests that during initial stages of decomposition losses of carbon and phosphorus are proportional, but that there are smaller relative losses of nitrogen due to interactions between soil organic matter and the physical soil environment (Tipping et al., 2016).Reduced losses of nitrogen during initial stages of decomposition would lead to elevated N : P during early stages of soil development compared to later stages.As the simplistic soil decomposition in ORCHIDEE (Parton et al., 1993) omits interactions between soil organic matter and the physical soil environment (Doetterl et al., 2015;Tipping et al., 2016), it fails to reproduce the strong influence of litter stoichiometry on the overall soil stoichiometry.
At the 4.1 Myr site, the C : N : P stoichiometry of the soil organic matter is at 158 : 12 : 1 relatively close to the observed C : N : P ratio of 215 : 10 : 1.The realistic phosphorus content of soil organic matter indicates that the relative contribution of biochemical mineralization is sufficiently well www.geosci-model-dev.net/10/3745/2017/Geosci.Model Dev., 10, 3745-3770, 2017 calibrated in the model.It has to be noted that recent findings indicate a preferential physical stabilization of nutrientrich soil organic matter (Tipping et al., 2016) and a role of phosphatases in rendering organic compounds available as a carbon source to microbes (Spohn and Kuzyakov, 2013).Both findings challenge the classical view of a primary control of biochemical mineralization on the soil organic matter phosphorus concentration (Walker and Syers, 1976).Thus the common calibration approach (Wang et al., 2010;Goll et al., 2012;Yang et al., 2014) might be shortsighted.
The nitrogen fixation rate of 2.25 ± 0.96 g m −2 at the 300year site lies within the range observed among sites during early (< 150-year) soil development in Hawaii (Crews et al., 2001) (Table 4).The lower fixation rates at the 4.1 Myr site are due to the high mineral nitrogen availability.As the regulation mechanisms of nitrogen fixation are elusive (Barron et al., 2009;Vitousek et al., 2013) we only account for direct product inhibition control, a common regulation mechanism in biological systems, via soil mineral nitrogen concen-(Eq.B2), omitting a direct influence of phosphorus availability on nitrogen fixation (Vitousek et al., 2013).
In summary, the model is able to capture biomass stocks and NPP well and -more importantly -captures the contrasting nutritional states of vegetation among the two sites, indicated by the foliage N : P ratio.As we prescribe a common parametrization of vegetation characteristics for both sites, the differences in the leaf stoichiometry are the emergent outcome of the process governing the access of plants to nutrients and their response to two contrasting situations of nutrient availability which originate solely from differences in phosphorus inputs, labile phosphorus sorption capacity, and organic matter accumulation rate.
The model fails to reproduce differences in the allocation of NPP to different tissue and tissue turnover between sites, due to insufficient plasticity in tissue turnover and biomass allocation.The simulations suggest that the recycling and storage of nutrients in ecosystems subject to periods of drought or other nutrient-unrelated declines in foliage and growth are an important source of nutrients for new growth.We show that the approach of calibrating biochemical mineralization rates using the soil organic matter stoichiometry is problematic, which is in line with growing evidence of the physical stabilization of soil organic matter (Doetterl et al., 2015;Tipping et al., 2016) and new insights into the functioning of phosphatase (Spohn and Kuzyakov, 2013).We further show that despite ORCHIDEE's intended application when designed, the model is able to capture to a large degree the general state of an ecosystem in an early stage of soil development.

Fertilization experiment
To evaluate the simulated response of vegetation to nutrient addition, we perform three nutrient treatment simulations per site: addition of either nitrogen (+N) or phosphorus (+P), or the combined addition of both (+NP).The annual addition rates of nitrogen and phosphorus of 10 g m −2 yr −1 are similar to the field experiments.
The model captures the signs of change in net primary productivity (NPP) to the fertilization treatments at both sites: at the 300-year site NPP strongly reacts to the addition of nitrogen, but not to the addition of phosphorus, while at the 4.1 Myr site NPP strongly reacts to the addition of phosphorus, but not to the addition of nitrogen (Fig. 4).For the 4.1 Myr site, the size of the simulated response ratio (defined as the NPP of the nutrient addition experiment divided by the NPP of the control experiment) is also comparable to the observed response ratio.At the 300-year site, the simulated positive effect of nitrogen deviates from the observation but is still within the range of uncertainty.However, the synergistic effects of combined addition are not captured by the model, as the model simulates a background phosphorus availability at the 300-year site which is high enough to support the nitrogen stimulated growth.
The 300-year site is accumulating organic material, in particular soil organic matter, and the accompanied immobilization of soil nutrients is the major driver of nutrient scarcity.This leads to extremely low mineral nitrogen concentration, whereas the high weathering release of phosphorus from minerals is sufficient to keep soil labile phosphorus concentration relatively high (Table 4).Therefore, vegetation reacts Geosci.Model Dev., 10,2017 www.geosci-model-dev.net/10/3745/2017/Table 5.The response of net primary productivity, leaf area index, nutrient use efficiencies, and foliar stoichiometry to nutrient addition at the 300-year old sites and the 4.1 Myr old site along the Hawaii chronosequence.Given are the relative difference between the nutrient addition experiments and the control simulation and, in parentheses, in the observed changes (Harrington et al., 2001;Ostertag, 2001) strongly to the addition of nitrogen at the young site.The lack of any stimulation of plant productivity in the model to phosphorus addition at the young site indicates an overestimation of plant available phosphorus, likely due to the omission of differences in the occlusion rate of soil labile phosphorus among sites, which tends to be much higher at the young site (Violette, unpublished data).At the 4.1 Myr site the remobilization of phosphorus from soil organic matter is the major source of phosphorus for vegetation, as the minerals are phosphorus depleted, leading to low soil labile phosphorus concentration (Table 4).Compared to the young site, a higher fraction of soil labile phosphorus is also adsorbed to soil particles and is thus not available to plants.Therefore, vegetation reacts strongly to the addition of phosphorus at the old site, but not to nitrogen addition.The increases in the NPP per leaf area (NPP/LA) are more pronounced than the increases in leaf area index when nitrogen stress is alleviated (300-year site), whereas the increases in NPP/LA are comparable to the increases in leaf area index when phosphorus stress is alleviated (4.1 Myr site) (Table 5), in both model and observations.Such a model behavior can be expected due to the lack of a direct effect of foliage phosphorus concentration on photosynthesis (Fig. 2).While the link between leaf nitrogen concentration and the carboxylation efficiency of photosynthesis (V cmax ) is well established (Kattge et al., 2009), the role of leaf phosphorus concentration in photosynthesis is less clear as nitrogen and phosphorus concentrations usually co-vary (Reich et al., 2009;Kattge et al., 2009;Domingues et al., 2010;Walker et al., 2014;Bahar et al., 2016;Norby et al., 2016).
The model fails to capture the drop in leaf area index when the non-limiting nutrient is added (which is phosphorus for the 300-year site and nitrogen for the 4.1 Myr site).The causes of the drop are unclear (Harrington et al., 2001), but could be related to an increase in grazing in the fertilized plots due to the higher nutrient content of foliage compared to the surroundings (Casotti and Bradley, 1991;Campo and Dirzo, 2003).
The model fails to capture the extent to which phosphorus use efficiency declines when phosphorus is added.This can be attributed to the omission of excessive plant uptake of nutrients (luxury consumption) which drives the observed reduction in the use efficiencies of the non-limiting nutrient (Harrington et al., 2001).As luxury consumption does not directly affect plant growth (Lawrence, 2001;Van Wijk et al., 2003) and is strongly species dependent (Lawrence, 2001), it is omitted in the model.
The model tends to overestimate increases in leaf nitrogen concentration in response to nitrogen addition, while increases in leaf phosphorus concentration in response to phosphorus addition are underestimated (Table 5).Biases in the availability of added nutrients to vegetation could be responsible for mismatch between observed and simulated responses.ORCHIDEE is prone to overestimating sorption losses when soil labile phosphorus concentration is substantially elevated above natural levels as the linear phosphorus sorption applied here cannot capture the increase in the dissolved fraction when the sorbed fraction approaches the maximum sorption capacity of the soil.The simulation setup of applying fertilizer evenly during the year, instead of in two fertilization events, underestimates nitrogen losses due to the nonlinear relationship between soil emissions and substrate availability (Shcherbak et al., 2014).However, the stoichiometric adjustments in ORCHIDEE are not process based and might be themselves the cause of the bias in the simulated response of leaf stoichiometry.
The sign of the responses of biochemical mineralization to nutrient addition mirrors the observed changes in potential phosphatase activity (Table 5).This indicates that the simplistic approach for biochemical mineralization applied in ORCHIDEE seems to capture the general behavior of phosphatase activity.However, it does not allow us to draw www.geosci-model-dev.net/10/3745/2017/Geosci.Model Dev., 10, 3745-3770, 2017 any conclusion about the overall significance of biochemical mineralization to total mineralization, which represents a major uncertainty for modeling phosphorus cycling (Goll et al., 2012;Yang et al., 2014;Reed et al., 2015).In summary, the model is able to capture the contrasting responses to the three fertilization treatments among the two sites.We find further that differences between sites in the underlying changes in LAI and NPP/LA are partly captured by the model: the alleviation of nitrogen stress at the 300-year site increases NPP/LA more strongly than LAI, while the alleviation of phosphorus stress increases NPP/LA and LAI equally (Table 5).Nonetheless, the model underestimates the nitrogen limitation of productivity at the 300-year site.The nitrogen capital of this site strongly depends on the nitrogen inputs and the efficiency at which nitrogen is retained in the ecosystem, which we could only roughly approximate.

Conclusions
Here, we present the implementation of a terrestrial phosphorus cycle and its interactions with the carbon and nitrogen cycle in the ORCHIDEE land surface model.The model accounts for effects of nutrient stress on tissue nutrient concentration, litter quality, root-to-shoot allocation, and photosynthesis.We further account for root phosphorus uptake and the movement of phosphorus in the soil volume as an additional constraint on soil phosphorus availability to plants, to reduce the sensitivity of plant phosphorus availability to the sorption dynamics which can only be poorly constrained from available data (Goll et al., 2012;Reed et al., 2015).
We evaluated the performance of the model at two sites of contrasting nutrient availabilities (Crews et al., 1995;Harrington et al., 2001;Ostertag, 2001).The model captures the different sensitivities of net primary productivity to nutrient addition among sites (Fig. 4).It further tends to reproduce differences in the leaf area index and leaf level productivity between the alleviation of nitrogen and phosphorus stress (Table 5).As we prescribed a common parametrization for all biological processes for both sites (Table 3), the contrasting response of vegetation to nutrient addition among sites is the emergent outcome of differences in the physico-chemical soil characteristics.
The model shows some deficiencies which can be linked to the lack of plasticity in the allocation of new growth to the different plant tissues and biomass turnover -a common issue in global models (De Kauwe et al., 2014).It further underestimates nutrient use efficiencies in general, primarily due to the overestimation of wood nutrient content and wood growth, but we cannot rule out that the nitrogen use efficiency of photosynthesis in ORCHIDEE which is derived from data of plants growing on a wide range of soil nutrient availability (Kattge et al., 2009) does not apply to the extreme environment found during early and late stages of soil formation.Data availability.Primary data and scripts used in the analysis and other supplementary information that may be useful in reproducing the author's work can be obtained by contacting the corresponding author.
Geosci.Model Dev., 10,2017 www.geosci-model-dev.net/10/3745/2017/We find that at both sites the underestimations of N prod of −59.4 and −44.7 % are related to the overestimation of tissue nitrogen content (see "observed C i " in Table A4).The effect of biases in simulated biomass on N prod differs between sites, as biomass stocks, and subsequently nitrogen stocks, are underestimated at the 300-year site, while they are overestimated at the 4.1 Myr site (Fig. 3).When biomass and nitrogen content are taken from observation, the bias in N prod is in general low due to the good agreement of simulated and observed NPP (Fig. 3).
The effects of biases in xc i and/or C i on P prod are comparable to their effects on N prod at both sites.
We calculated the residence time of nutrients (X MRT,i ) for classes of tissues sharing similar stoichiometry in OR-CHIDEE, namely leaf, coarse root and stems, and fine roots: where NPP i is the fraction of NPP being allocated to tissue class i.The equation can be simplified to showing that the residence time X MRT,i is equal to the residence time of carbon.Table A5 shows the residence times of the different tissue classes as well as C MRT , N MRT , and P MRT calculated as described in the method section of the main paper.
The observed residence times of carbon, nitrogen, and phosphorus (Table A5) are longer at the 4.1 Myr site than they are at the 300-year site, whereas the model simulates an opposite pattern.As the constant tissue turnover rates are used in ORCHIDEE, the model is not able to reproduce differences among sites.The slightly longer residence times at the young site can be attributed to the transient state of vegetation in which it is still accumulating biomass, while the biomass at the old site reached a stable state.The model consistently underestimates the residence time of phosphorus P MRT , which can be attributed to an underestimation of the phosphorus content of stems and coarse roots (not shown).

Figure 1 .
Figure1.Pools and fluxes of the phosphorus component of ORCHIDEE: phosphorus enters ecosystems in the form of dissolved labile phosphorus from chemical weathering of minerals and atmospheric deposition.Dissolved labile phosphorus is the sole source for plants and microbes and can be reversibly adsorbed onto soil particles or lost by leaching.Adsorbed labile phosphorus can be transformed into more recalcitrant forms which are irreversibly lost to biota.When plants take up phosphorus from soils it enters the plant labile phosphorus pool, from which it is allocated to growing plant tissues.The plant labile phosphorus concentration is buffered by a reserve pool which serves as a long-term storage to buffer seasonal variations in phosphorus demand and supply.When plant tissue is shed, part of the phosphorus is recycled (broken lines), while the rest enters the litter pools, from where it is either transformed into soil organic matter or mineralized.

Figure 3 .
Figure 3.Comparison between simulations (mean ± SD) and observations (mean ± SD -if available) at the 300-year old site (a) and the 4.1 Myr old site (b) along the Hawaii chronosequence.Net primary productivity (NPP) is given in g (C) m −2 yr −1 ; the standing stocks of foliage, stems (including coarse roots), fine roots, and soil organic matter (SOM) are given in g (C) m −2 .The SD in the case of simulations indicates the inter-annual variability, not model error.

Figure 4 .
Figure 4. Comparison of simulated and observed responses of net primary production to fertilization at the 300-year old site (a) and the 4.1 Myr old site (b) along the Hawaii chronosequence.The response ratio is the measured or modeled plant production in the fertilizer treatment divided by its value under unfertilized conditions.The bars represent the measurement uncertainty and the annual variability in simulations, respectively.
Code availability.The ORCHIDEE model version used here is a development branch of ORCHIDEE, which is open source.The SVN version of the code branch is https://forge.ipsl.jussieu.fr/orchidee/browser/branches/ORCHIDEE-CN-P revision 4520 from 21 July 2017.Please contact the corresponding author for the code of ORCHIDEE-CN-P if you plan an application of the model and envisage longer-term scientific collaboration.

Table 2 .
Parameters of the model.

Table 3 .
The site conditions prescribed for the simulation for the 300-year old site and the 4.1 Myr old site along the Hawaii chronosequence.
. We used a Student's t test to test whether nutrient addition resulted in a significant difference in the respective variables.If significant (p > 0.05), the values are in bold.

Table A2 .
Parameters of the dynamical weathering routine.

Table A3 .
Additional variables of the dynamical weathering routine.

Table A4 .
The bias (simulated-observed) in carbon production per biomass nutrient (N prod , P prod ).From simulations (default) and as diagnosed by substituting simulated with observed nutrient content (observed xc i ) and/or biomass stocks (observed C i ).

Table A5 .
The overall residence times of carbon (C MRT ), nitrogen (N MRT ), and phosphorus (P MRT ) in biomass (excluding storage), as well as the approximated residence time of all elements (C, N, P) in foliage (X MRT,leaf ), coarse roots, and stems (X MRT,wood ), as well as fine roots (X MRT,root ).The latter three tissue classes consist of tissue with similar stoichiometry and thus the residence time of all elements is the same.