A global model of carbon, nitrogen and phosphorus cycles for the terrestrial biosphere

. Carbon storage by many terrestrial ecosystems can be limited by nutrients, predominantly nitrogen (N) and phosphorus (P), in addition to other environmental constraints, water, light and temperature. However the spatial distribution and the extent of both N and P limitation at the global scale have not been quantiﬁed. Here we have developed a global model of carbon (C), nitrogen (N) and phosphorus (P) cycles for the terrestrial biosphere. Model estimates of steady state C and N pool sizes and major ﬂuxes be-tween plant, litter and soil pools, under present climate conditions, agree well with various independent estimates. The total amount of C in the terrestrial biosphere is 2767 Gt C, and the C fractions in plant, litter and soil organic matter are 19%, 4% and 77%. The total amount of N is 135 Gt N, with about 94% stored in the soil, 5% in the plant live biomass, and 1% in litter. We found that the estimates of total soil P and its partitioning into different pools in soil are quite sensitive to biochemical P mineralization. The total amount of P (plant biomass, litter and soil) excluding occluded P in soil is 17 Gt P in the terrestrial biosphere, 33% of which is stored in the soil organic matter if biochemical P mineralization is modelled, or 31 Gt P with 67


Introduction
Simulations using global climate models with a fully coupled carbon cycle showed that warming Could reduce the net carbon storage in the terrestrial biosphere globally, resulting in an increase in atmospheric CO 2 concentration and further warming of 0.1 to 1.5 • C by 2100 (Friedlingstein et al., 2006).However there are considerable uncertainties in those predictions.For example, none of those models explicitly included nutrient limitations and their responses to climate and higher (CO 2 ).Both field measurements and theoretical studies have shown that nitrogen limitation can have a significant influence on how the carbon cycle will respond to increasing (CO 2 ) (Luo et al., 2004) and warming (Medlyn et al., 2000).This is also supported by recent studies (Sokolov et al., 2008;Churkina et al., 2009;Thornton et al., 2009;Wang and Houlton, 2009;Zaehle et al., 2010).
Globally N and P are the most common nutrients limiting plant growth and soil carbon storage (Vitousek and Howarth, 1991;Aerts and Chapin, 2000).A number of global biogeochemical models have been developed to account for N limitation on the productivity of and C uptake by the terrestrial biosphere (Parton et al., 1987;McGuire et al., 1995;Thornton et al., 2009;Xu-Ri and Prentice, 2008;Zaehle et al., 2010;Gerber et al., 2010), but only the CENTURY model (Parton et al., 1987) simulates biogeochemical cycles of C, N and phosphorus (P) and its P cycle submodel has yet to be applied globally.There are some strong reasons why the P cycle should be included in global models for studying the interactions between climate and biogeochemical cycles: (1) both theory and experiments suggest that much tropical forest and savannah are phosphorus limited (Aerts and Chapin, 2000), and tropical forests and savannahs account for about 40% of global vegetation biomass (Saugier et al., 2001) and 45% of global terrestrial net primary productivity (Field et al., 1998) ; (2) a recent study by Houlton et al. (2008) showed that biological N fixation, the largest N input to the un-managed terrestrial ecosystems at present is closely related to phosphatase production in the tropics; (3) responses of N and P cycles to climate, increasing atmospheric (CO 2 ) and human activities can be quite different because of the different biogeochemical controls on N and P cycles in the terrestrial biosphere (Vitousek et al., 1997).For example, the external input to the unmanaged ecosystems is dominated by N fixation for N, but by weathering and dust deposition for P for most unmanaged lands.Loss from the unmanaged ecosystems is dominated by gaseous fluxes via denitrification or leaching for N and by phosphate leaching for P. Misrepresenting nutrient limitation in the tropics may lead to incorrect predictions under future climate conditions.An early study showed that the relative response of leaf photosynthesis to elevated (CO 2 ) is smaller when plant growth is P limited (Conroy et al., 1990) as compared to the response under N-limited conditions and (4) some terrestrial ecosystems may shift from N limitation to P limitation under high N input (Perring et al., 2008) or future climate and higher (CO 2 ) conditions (Menge and Field, 2007;Matear et al., 2010).
The objectives of this study are (1) to develop a global biogeochemical model of C, N and P cycles for the terrestrial biosphere for use in a global climate model or earth system model; (2) to construct steady state C, N and P budgets for the terrestrial biosphere for the 1990's using available information of plant biomass, litter fall rate and soil C and N and estimates of P for different soil orders; and (3) to provide a quantitative estimate of the extent of N and P limitations and their uncertainties on plant productivity globally under the present conditions.Our model calibration strategy assumes that all fluxes are in steady state in the 1990s and the limitation of this strategy is discussed later (see Sect. 6).
While a reasonable amount of information is available for the pools and fluxes of the terrestrial carbon and nitrogen cycles (Post et al., 1982, 1985, Field et al., 1998 for example), global datasets are scarce for the P cycle and highly uncertain.Spatially explicit estimates of soil P amount are not yet available globally.Even the estimates of global total amount of soil P vary widely from, for example 200 Gt P (Jahnke, 1992) to 40-50 Gt P (Smil, 2000), due to different assumed mean P content (0.1 or 0.05%) and soil thickness (60 or 50 cm).To overcome the data limitation, we will use a well calibrated carbon cycle model (CASA') as our carbon cy-cle submodel and modeling framework and make use of observed coupling among all three cycles to constrain the nutrient (N and P) pools and fluxes.
The model CASA' was developed from the CASA model (Randerson et al., 1997) and has been used in studying the carbon-climate feedback globally (Fung et al., 2005).CASA' uses NPP from a coupled land surface model whereas CASA is an offline model using satellite derived NPP (Randerson et al., 2009).We have added N and P cycles to CASA' by adapting the N and P cycle model developed by Wang et al. (2007) and Houlton et al. (2008) from single litter and soil organic matter (SOM) pools to the multiple litter and SOM pools used by CASA'.This is important because a recent study showed that multiple-pool representation is required for studying the response of soil respiration at decadal or century time scale (Knorr et al., 2005).
Because the carbon cycle of CASA' model has been well calibrated (Randerson et al., 1997) and has been used in several previous studies (Fung et al., 1997;Randerson et al., 2002 for example), we address the question here of what sizes the nutrient pools and fluxes should be for the global C cycle as represented in the CASA' model for the 1990's.To estimate the pool sizes and fluxes of N and P, we drive the model using the spatially explicit estimates of monthly nutrient-unlimited NPP for the 1990's as input to our model, and calculate nutrient-limited NPP, and nutrient limitation factors relative to that in the 1990's for each land point at steady state.The steady state assumption is used to reduce the dependence of the estimates of N and P pools and fluxes on their initial estimates for which we currently have little global-scale spatially explicit information.The couplings of the three cycles, as represented by our model, are calibrated using independent estimates globally, such as leaf N:P ratio and fraction of P in different soil pools.The modelled pools and fluxes are then compared with estimates from other studies.
In Sects.2-4, we describe the model, model calibration, and model evaluation under present climate conditions against independent estimates of various pool sizes and biogeochemical fluxes at global scales.Section 5 describes the predicted nutrient limitation globally under the present conditions.Section 6 discusses limitations of the present study, and future studies to address those limitations.

Model description
The pools used to represent the C, N and P cycling through the terrestrial ecosystem in plants, litter and soil are shown in Fig. 1.Plants are divided into leaf, wood and root pools, litter into metabolic litter, structural litter and coarse woody debris pools and soil into microbial biomass, slow and passive pools.The turnover rate depends on soil temperature, moisture and texture for litter and soil pools (Randerson et al., 1997) or biome for plant pools.There is one additional 59 deposition (N and P), weathering (P), fixation (N) and fertilizer addition (N and P), output in red is loss by leaching or gaseous loss from the ecosystem.Plants take N from the inorganic N pool and P from the labile P pool in soil.), soil (yellow brown) into microbial biomass, slow pool and passive pool.One inorganic soil mineral N pool and three other P pools are also represented.Arrows between the pools represent the direction of C, N and P flow between pools.For N and P, external inputs are deposition (N and P), weathering (P), fixation (N) and fertilizer addition (N and P), output is loss by leaching or gaseous loss from the ecosystem.Plants take N from the inorganic N pool and P from the labile P pool in soil.
pool for N (inorganic N (NO − 3 + NH + 4 ) in the soil) and three additional P pools (labile, sorbed and strongly sorbed P) in our model.Change in a pool size with time is governed by a differential equation that is numerically integrated daily.We shall present an overview of each of the three cycles and their interactions in the following sections.A detailed description including key equations and parameter values is given in the appendices.A full list of symbols and their definition are provided in Appendix A.

Carbon cycle
The carbon cycle is based on CASA' model (Fung et al., 2005).We reduced the number of carbon pools by combining surface litter with soil litter, and surface microbial biomass with soil microbial biomass.This gives three discrete pools in the litter: structural, metabolic and coarse woody debris pools and three organic pools in the soil: microbial biomass, slow and passive pools.The fluxes between different pools are modeled as in CASA'.Details are given in Appendix B.
Transfer coefficients from plant pool i to litter pool j , b j,i and from litter pool j to soil pool k, c k,j are calculated as in CASA' model (Fung et al., 2005).Turnover rates of litter carbon (µ j ) or soil carbon (µ k ) are a function of substrate quality (lignin:N ratio), soil temperature, moisture and soil texture (Randerson et al., 1997).The turnover rate of leaves is calculated as a function of leaf age (Arora and Boer, 2005), and the turnover rates of woody tissue or fine roots are constant for each biome, but vary with biome type (see Table 1).
Because the N:C ratios of litter pools are much lower than those of soil, decomposition of litter carbon can be limited by available soil mineral N. When litter decomposition is not N limited, decomposition of litter and soil is limited by the amount of substrate, not its quality.When litter or soil carbon is decomposed, some of the decomposed carbon is respired as CO 2 .Heterotrophic soil respiration is calculated as the sum of the respired CO 2 from the decomposition of all litter and soil organic C pools.We assumed that the storage change of gaseous CO 2 in the soil is negligible, therefore the surface CO 2 flux is equal to CO 2 production in the soil.The difference between NPP and soil respiration is net ecosystem C exchange (NEE) between the land surface and atmosphere.
Input to the carbon cycle includes nutrient unlimited NPP, and initial carbon pool size, output are nutrient limited NPP, soil respiration, NEE and model pool sizes.Nutrient unlimited NPP can also be provided by a global land surface model when CASACNP is coupled to a global climate model.In this study, we used the scaled NPP from CASA simulation as the nutrient unlimited NPP (see Sect. 2.5).

Nitrogen cycle
The nitrogen cycle is based on the model developed by Parton et al. (1987) and Wang et al. (2007).Similar to the C cycle, the change in N in each pool is governed by a differential equation (see Appendix C).An additional mineral N pool in soil is also represented as only mineral N is assumed to be taken up by plants.We do not include uptake of organic N in soil by roots (Schimel and Bennett, 2004).Ammonia volatilization is not modeled as it usually occurs when soil pH is above 8 (Freney et al., 1983), and the fraction of land with pH>8 is very small globally (Batjes, 1996).
We do not explicitly model the processes of nitrification and denitrification.Therefore our model will need further improvement in the future.In our model gaseous N loss is assumed to be proportional to net N mineralization based on the "holes-in-the-pipe" idea (Firestone and Davidson, 1989) and the rate of leaching loss is proportional to the soil inorganic N pool size.Leaching loss of soil organic matter is not included in our model.
The nitrogen cycle is closely coupled to the carbon cycle; carbon decomposition and gross N mineralization is coupled by the N:C ratios of the substrates (compare equations B2 and B3 with equations C2 and C5).Net N mineralization rate (F n,net ) is the difference between gross N mineralization (F n,gr ) and N immobilization (F n,im ).When net mineralization rate is negative (gross N mineralization < N immobilization), and the additional amount of mineral N required by N immobilization can not be met by the amount of mineral N available, the litter carbon decomposition rate is reduced (see Appendix C Eq. C12 for m n ).
Nitrogen uptake by plants is modeled as a function of soil mineral N pool size, and the demand by plant growth (Eq.C7), similar to the TEM model (Melillo et al., 1993).The nitrogen demand is a product of maximal N:C ratio and NPP allocated to each plant pool minus the amount of  Wetland, urban land and land ice in the IGBP biome classification are not included in our simulations.The mean N:C ratio of leaves is based on estimates from the Glopnet datasets for each biome (Wright et al., 2004) and mean P:C ratio for each biome is calculated from the mean leaf N:C ratio from Glopnet dataset and the estimated N:P ratio from this study for each biome.Minimal and maximal leaf N:C or P:C ratios are assumed to be 0.8 and 1.2 times the mean leaf N:C or P:C ratios for each biome.NPP allocation coefficients during steady leaf growth (a leaf , a wood , a root ) and mean residence time of plant tissue (1/µ i ) are based on the CASA model.Estimates of N:C and P:C ratios are based on Weedon et al. (2009) for woody tissues and on Gordon and Jackson (2000) for roots.Parameters leaf N:P ratio, x npmax , v pmax and soil C:N are estimated for each biome during model calibration.The numbers in brackets for leaf N:P ratio are the one standard error of the mean from our calibration.Leaf N:P ratio is not fixed in model simulations.resorbed N from that pool.When the uptake is greater than the minimal demand, the amount of uptake nitrogen allocated to each pool is in proportion to the demand.During senescence, some fraction of plant tissue nitrogen is resorbed to live tissue, and the remaining goes to the litter pool.Leaf and root litter are partitioned into metabolic litter and structural litter.The N:C ratio is fixed for structural litter (=1/125) but variable for metabolic litter.Woody litter goes to the coarse woody debris pool directly.Only the N:C ratio of soil organic matter and structural litter pools are fixed.N:C ratios of all plant pools are allowed to vary within prescribed ranges (see Table 1).
Input of N to the model includes atmospheric N deposition (both wet and dry), N fertilizer application, N fixation (both symbiotic and asymbiotic) and output includes N leaching and gaseous loss.

The phosphorus cycle
The phosphorus cycle is based on the model of Wang et al. (2007) and Houlton et al. (2008).The differential equations used to describe the rate of change of each pool are presented in Appendix D. Three differential equations are used to represent the dynamics of labile, sorbed and strongly sorbed phosphorus in soil.The P:C ratio for the three different plant pools.can vary within a given range for each biome, therefore N:P ratios of plant and litter pools can vary.
The N:P ratios of the newly formed soil organic pools are fixed.However the N:P ratios of the slow and passive pools will change as P in these two pools can be mineralized both biologically and biochemically.The biological P mineralization is the same pathway as N mineralization by microbial activities, and the rate of gross biological P mineralization is calculated as the carbon decomposition rate divided by the P:C ratio of the substrate.P immobilization rate is calculated as the N immobilization rate divided by the N:P ratio of different soil pools.The N:P ratio of the newly formed soil organic pool is 4 g N (g P) −1 for microbial biomass (Cleveland and Liptzin, 2007) and 7 g N (g P) −1 for the slow and passive pools for highly weathered soil orders and 5 g N (g P) −1 for the other soil orders (Crews et al., 1995).
Phosphorus in the slow and passive soil pools can also be mineralized biochemically (McGill and Cole, 1981).Therefore the N:P ratios of the slow and passive pools will increase until a steady state is reached when the P fluxes into those two pools through P immobilization (biologically only) are equal to the rates of P being mineralized (both biologically and biochemically) from those pools.Biochemical mineralization is modeled as a function of soil organic P, the N costs of P uptake and phosphatase production, and maximal specific biochemical P mineralization rate (see Wang et al., 2007;Houlton et al., 2008).
We do not model the biochemical P mineralization of litter P, as turnover rates of the litter pool are much faster than those of the slow and passive soil pools, and all P in the litter will be mineralized biologically if they are not mineralized biochemically within a few years.We do not distinguish the phosphatase production by roots from that by soil microbes in our model.
We assumed that the labile P pool is equilibrated with the sorbed P within days.The relationship between the amount of labile P and sorbed P is described using the Langmuir equation (Barrow, 1978;Lloyd et al., 2001;Wang et al., 2007).Inputs to the labile P pool are net biological P mineralization and biochemical P mineralization, P weathering, dust deposition and P fertilizer addition.Only labile P can be taken up by plants.
Some of the sorbed P can enter the strongly sorbed P pool that is not exchanged readily with the labile P; the rate of sorbed P to strongly sorbed P is assumed to be proportional to the amount of sorbed P in the soil.The flux from the strongly sorbed P pool to occluded P pool that is not available to plant or soil microbes at a time scale of decades to a century is not represented in our model.Including the dynamics of occluded P pool will significantly increase the computation with little impact on the simulated processes we are interested in here at decade or century scales.
Because of the biochemical P mineralization, the P cycle in the soil can become quite decoupled from C and N cycles in the soil (McGill and Cole 1981).However a recent study by Houlton et al. (2008) showed that the N cycle may be significantly coupled to the P cycle in some tropical soils, as the N fixation is dependent on the rate of biochemical P mineralization and N and P cycles in the N-limited tropical soils can be strongly coupled.Since we do not simulate N fixation explicitly, this coupling between N and P cycle has not yet been included in our present model.
Inputs of P to the ecosystem are weathering, deposition and fertilizer application.Outputs are the leaching loss of labile P and loss of strongly sorbed P to the occluded P.

Nutrient limitation on net primary productivity
We model NPP as a function of two nutrient limitation factors.That is where F cmax is the nutrient unlimited NPP (g C m −2 day −1 ) and x npleaf is the nutrient concentration limiting factor, and x npup is the nutrient uptake limiting factor.x npleaf is calculated as and where n leaf and p leaf are the N:C (g N/g C) and P:C (g P/g C) ratio of the leaf biomass, k n and k p are two empirical constants.
Nutrient uptake limiting factor, x npup is calculated as and x p,up = min 1, where N min is the amount of mineral N in soil (g N m −2 ), and P lab is the amount of labile P in soil (g P m −2 ), t is the time step of model integration (=1 d), F n,upmin and F p,upmin are the amount of minimal N and P uptake required to sustain a given NPP.Therefore nutrient uptake limiting factor will become less than 1 when the available nutrients (N or P) amount is less than the minimal amount of nutrient required by plants for a given NPP.Equation (2) states that NPP is limited by N or P. When NPP is N limited (x n,leaf < x p,leaf ), increasing p leaf will not reduce the N limitation, and vice versa.This is supported by the results from fertilizing experiments (Vitousek, 2004).Equation ( 3) is used because both photosynthesis and plant respiration increase with n leaf , and the increase in photosynthetic rate per unit leaf N is slower than that of respiration per unit leaf N at higher leaf N (Kattge et al., 2009, Reich et al., 2005).A similar model is commonly used in estimating the response of NPP to nitrogen limitation (Melillo et al., 1993;McMurtrie, 1991).
Very few measurements are available on the responses of photosynthesis or respiration to p leaf .Some earlier measurements, summarized by Lloyd et al. (2001), suggest that the response curve of leaf photosynthesis to p leaf has a similar shape to that for n leaf , which is also consistent with a more recent study on tropical grasses (Ghannoum et al. 2008).A study by Kattge et al. (2009) also found that the estimated maximum carboxylation rate per unit leaf nitrogen of the leaves of tropical forests on the phosphorus-poor oxisol soils is lower than those on other soils.Meir et al. (2000) found that the respiration of leaves of tropical trees was better correlated with p leaf than with n leaf .Observations from longterm fertilization experiments also show that both N-limited and P-limited forests responded to applications of the limiting fertilizer by increasing canopy leaf area, radiation use efficiency and foliar nutrient concentration (Harrington et al., 2001).For these reasons, we used the same function with different model parameters for estimating x n,leaf and x p,leaf .
In this study, we assume that k n =0.01 g N (g C) −1 , based on the results of Linder and Rook (1984) and k p =0.0006 g P (g C) −1 .It has been suggested that NPP is N limited when leaf N:P (on mass basis) <14 and is P limited when leaf Y. P. Wang et al.: A global model of carbon, nitrogen and phosphorus cycles N:P>16 based on broad-scale geographic variations of leaf N:P ratios (Koerselman and Mueleman, 1996).The value of k p is chosen so that NPP is limited by N (x n,leaf < x p,leaf ) when n leaf /p leaf <16 (g N/g P) and otherwise NPP is limited by P. As discussed in Aerts and Chapin (2000), this is a first approximation for studying nutrient limitation at broad scales.Globally the variation of leaf N:P ratio is found to be consistent with the expected nutrient limitation on NPP (Reich and Oleksyn, 2004).However species composition and other factors also likely affect the nutrient limitation within an ecosystem (e.g.Townsend et al., 2007).
Our model, as described so far, calculates nutrient-limited NPP by accounting for nutrient feedback on NPP when leaf N:C or P:C changes or soil nutrient supply cannot meet plant demand.One of the objectives of this study is to estimate the nutrient limiting factors under present conditions.Our approach here is to estimate the monthly nutrient unlimited NPP, F cmax , by multiplying the monthly NPP estimates of Randerson et al. (1997) at 1 • by 1 • spatial resolution for the 1990's by a biome-specific parameter x npmax (see 2.5).We run the model to steady state using the monthly F cmax at 1 • by 1 • as an input.At steady state, the nutrient limiting factor is equal to x npleaf , because x npup =1.However, when the external environment is changed, such as through an increase in atmospheric (CO 2 ), the nutrient uptake may limit NPP, and progressive nutrient limitation can occur (Luo et al., 2004).
Because we allow both N:C and P:C ratios of leaf biomass to vary within their prescribed ranges for each biome (see Table 1), the modelled leaf N:P ratio, ecosystem NPP will vary depending on the nutrient supply and demand.

Values of model parameters
The model has a total of 31 pools: 9 C pools, 10 N pools and 12 P pools.The C cycle of the model was calibrated using global data of (CO 2 ), 14 CO 2 (Randerson et al., 1997(Randerson et al., , 2002) ) and used for global studies (Fung et al., 1997(Fung et al., , 2005)).We used the same turnover rates and transfer coefficients for all litter and soil pools as Randerson et al. (1996).
For each biome, we prescribed the ranges (minimal and maximal) of C:N ratios of leaves based on data compiled in Glopnet (Wright et al., 2004), the ranges of C:N:P ratio of wood based on the results of Weedon et al. (2009), and the ranges of C:N:P ratio of roots based on Gordon and Jackson (2000).We calculated the range of C:P ratio of leaves from the C:N ratio and the estimated N:P ratio of leaves for each biome (see Table 1).We estimated the leaf N:P ratio for each biome from calibration using the empirical relationships between leaf N:P ratio and latitude by Hedin (2004).The actual C:N:P ratios of all plant pools during model integration will vary from point to point, depending on the available soil nutrients (N and P) for plant uptake.In this study we allow the C:N and C:P ratios of all plant pools to vary within their prescribed ranges.When the minimal nutrient demand can not be met by the available nutrients in soil, nutrient uptake limi-tation will reduce the nutrient-unlimited NPP, therefore both N:C and P:C ratios of all plant pools will not fall below their respective prescribed minima.
C:N:P ratio of the structural litter pool was fixed at 3750 (g C):25 (g N):1 (gP) for all biomes, and the C:N:P ratio of the metabolic litter pool was allowed to vary, depending on the quality of litter input.C:N ratios of soil organic pools are fixed for each biome.C:N:P of soil microbial biomass were fixed at 32 (g C):4(gN):1 (gP) for all biomes based on the estimate of Cleveland and Liptzin (2007).C:N ratio of slow and passive soil organic matter was estimated for each biome by calibrating the soil N estimate for each 2.5 • latitudinal bands against the estimate of Post et al. (1985).N:P ratios of newly formed slow and passive soil organic matter are assumed to vary with soil order, being 7 for the highly weathered soils and 5 for the rest (see Table 2) based on results of Crews et al. (1995).As soil ages, the N:P ratio of slow and passive pools will vary, depending on soil P biochemical mineralization.
To use Eq. ( 1) in our simulations, the value of x npmax is required for each biome and is estimated as follows.Using the prescribed monthly NPP of Randerson et al. (1997) as F c (F c =F c,1990 ) in our model, we ran the model to steady state to determine the pool sizes.Using these estimates of all pool sizes at steady state as the initial pool sizes, we ran the model again with F c =F c,1990 and calculated x npmax as where N c is the number of cells for each biome and A is cell area.That is equivalent to assuming that the mean biome x np is equal to 1 under steady state in 1990's.The value of x np can be considered as the nutrient limitation relative to the present conditions.Values greater than 1 indicate that the nutrient limitation is less than that under present conditions, and vice versa.Values of x npmax are listed in Table 1 for each biome.Because of variation of leaf N:C and P:C ratios within a biome, the nutrient limitation (x np ) can differ from 1 for some grid cells.Two parameters, k plab and s pmax affect the partitioning between labile P and sorbed P at equilibrium, and vary with soil order.Based on the estimates of different fractions of labile P, sorbed P and strongly sorbed P for different soil orders by Cross and Schlesinger (1995), we tuned these two parameters for each soil order using a nonlinear parameter estimation technique (Wang et al., 2009).Because of the strong Correlation between s pmax and k plab , we restricted the value of s pmax to be between 50 and 100% of the total inorganic soil P in the optimization.The value of k plab is consequently dependent on the range of s pmax as well as the values of the other parameters.
The biochemical P mineralization rate, F p,tase affects the model estimate of the fraction of organic P in soil, and is modelled as a function of the maximal specific biochemical P mineralization rate (v pmax ), the N cost of P uptake (λ pup ) and N cost of phosphatase production (λ ptase =15 gN/gP) (Treseder and Vitousek, 2001).Parameter v pmax is tuned to match the fraction of organic P in soil for each soil order.λ pup is 25 g N (g P) −1 for tropical evergreen broadleaf forests and savannahs, and is equal to 40 g N (g P) −1 for all other biomes, based on the simulation results of Houlton et al. (2008).For values of other parameters, see Tables 1 and  2.

Datasets
Three different kinds of datasets are used in this study: input data, data for model calibration and data for model evaluation.

Input dataset
We calculated the nutrient unlimited NPP for each grid cell by multiplying the monthly NPP from Randerson et al. (1997) for the 1990's by a biome-specific constant (x npmax ), and use them as input to the model.N input includes deposition, fertilizer application, fixation.We used the spatially explicit estimates of N deposition for 1990's by Dentener (2006) and N fixation by Wang and Houlton (2009) for the present climate conditions for N input.Global N input is 0.142 Gt N year −1 from fixation and 0.069 Gt N year −1 from deposition in 1990's, both are spatially explicit at 2 • by 2 • globally.Global N input from fertilizer application is taken as 0.086 Gt N year −1 (Galloway et al., 2004) and is distributed uniformly within the cropland biome.P inputs include fertilizer application, dust deposition and weathering.Global P input from fertilizer application is 0.014 Gt P year −1 in the 1990's (Smil, 2000) and is also distributed uniformly within the cropland biome.Spatially explicit P input from dust deposition is from the model output by Mahowald et al. (2008), and is 0.0007 Gt P year −1 globally in the 1990's.No spatially explicit estimates of global P weathering rates are available, and the few estimates of P weathering from different sites are highly variable (Newman, 1995).Estimates of P weathering rates from sites along a soil-age gradient have been used to estimate soil phosphorus content for relatively wet regions globally (Porder and Hilley, 2010).In this study, we estimated P weathering rate by dividing all 12 soil orders (Fig. 2) into four groups according to their weathering status.Soil orders within each weathering status were assigned a constant P weathering rate (Table 2).Based on the range of P weathering rates from soils along an age gradient in Hawaii (Chadwick et al., 1999), we assign the rate of 0.05 g P m −2 year −1 to the least weathered soils, such as Entisol, and low values to  1) and USDA soil order map (lower panel).Numbers 1 to 12 in the lower panel correspond to soil orders of Alfisol, Andisol, Aridisol, Entisol, Gellisol, Histosol, Inceptisol, Mollisol, Oxisol, Spodosol, Utisol and Vertisol, respectively.most weathered soils, such as Utisol (0.005 g P m −2 year −1 ) and Oxisol (0.003 g P m −2 year −1 ).These values are consistent with the estimated P weathering rates varying from 0.005 to 0.05 g P m −2 year −1 by Newman (1995) for different soils worldwide, and 0.007 g P m −2 year −1 for the highly weathered soils in the Amazon basin by Gardner (1990) for the residual soils (saprolites).The P weathering rate for the intermediate weathered soil orders was varied to give a global total P weathering rate of about 0.002 Gt P year −1 (Filippelli, 2002).

Datasets for model calibration
We did not carry out a rigorous calibration of all model parameters because insufficient global data are available, particularly for the P cycle.However components of the model that are important for nutrient limitation are calibrated using the estimates of soil N (Post et al., 1985), leaf N:P ratio (Hedin, 2004) and the estimates of the fractions of soil P in different pools (Cross and Schlesinger, 1995).The range of the leaf C:N ratios (dn leaf ) (prescribed for each biome) and the uncertainties of leaf P (dp leaf ) that are calculated as the product of leaf N:C ratio and P:N ratio (from calibration) are

Datasets for model evaluation
We used a number of datasets for evaluating the modeled pool sizes and fluxes.These datasets are: global vegetation biomass data (Olson et al., 1985), soil carbon pool size (Post et al., 1982), estimates of litter production (Matthews, 1997), global leaching and gaseous N losses, P leaching (Seitzinger et al., 2006) for N and P fluxes.These datasets are chosen because they are derived either directly from field observations or based on empirical relationships that are estimated from the field observations.Some of the datasets, such as vegetation biomass (Olson et al., 1985) and litter C production (Matthews, 1997) have spatially explicit information.However, as the biome classifications used by those authors are different from the IGBP biome classification we used in this study, spatially explicit comparisons could be misleading.Instead, we aggregated the spatially explicit estimates by the IGBP biome type, or by latitude for comparing with our estimates.Outputs from some other process-based models are also compared with our estimates.

Model integration
The model integration time step is one day.Meteorological inputs required for the model include daily surface air temperature, soil temperature and moisture.The daily meteorological forcing was generated using the CSIRO Conformal Cubic Atmosphere Model, CCAM, (McGregor and Dix, 2008) with the CSIRO Atmosphere and Biosphere Land Exchange (CABLE) land surface scheme (Wang and Leuning, 1998;Kowalczyk et al., 2006) at a spatial resolution of approximately 220 km globally.CCAM was run using sixhourly NCEP reanalysis for 1990 to 1997 (Kalnay et al., 1996) to produce daily mean air temperature, soil temperature and soil moisture in the rooting zone.
By reusing the daily forcings from 1990 to 1997, we ran the model to steady state.Steady state is considered to have been reached when the relative changes in total pool sizes of C, N or P per land point are less than 0.001% per year.All results reported here are for steady state in the 1990's only.Mass balances of all three cycles are achieved at every time step during the model integration.

Calibration of the biome -specific leaf N:P ratio, soil C:N ratio, and soil phosphorus fractions
Three datasets are used to calibrate leaf N:P ratio, C:N ratio of slow and passive soil pools, and the three soil parameters that affect the partitioning of soil P into different pools.
We calibrated the leaf N:P ratios of different biomes using the empirical relationship between leaf N:P ratio and latitude derived from field observations by Hedin (2004) (Fig. 3).Mean leaf N:P ratio by the calibrated model is above 16 g N/gP within the tropical region (about 15 • north or south of the equator), and less than 14 g N/g P in the region about 30 • away from the equator (see Table 1).The estimates of leaf N:P ratio by two other studies (Reich and Oleksyn, 2004;Kerkhoff et al., 2005) also fall within the range of the variation of leaf N:P ratio with latitude by our calibrated model.We then used the leaf N:P ratio as estimated here to calculate the range of C:P ratios for each plant biomass pools from the estimated C:N ratios from Glopnet dataset (see Table 1).In all CASACNP simulations presented here, we allow C:N and C:P ratios of each plant biomass pool to vary within their prescribed ranges, so N:P ratio of each plant biomass pool is not fixed.

61
(yellow green).The error bars represents the one standard error of the mean leaf N:P estimate by CASACNP within each latitudinal band.The biome mean leaf N:P ratio was calibrated against the relationship of Hedin (2004).Assuming that the C:N ratios of slow and passive soil organic matter are the same, we estimated that ratio for each biome by minimizing the squared difference between the C:N ratio of the modelled soil organic pools and the C:N ratio calculated from the estimates of soil C and N pools of Post et al. (1982Post et al. ( , 1985) ) for each 2.5 • latitudinal band (Fig. 4).Results are presented in Table 1.The modelled soil C and N pool sizes for the high northern latitudes (65 • N to 75 • N) are much smaller than the estimates by Post et al. (1982Post et al. ( , 1985)), because our model does not include wetland which has a very high content of soil organic matter (see Post et al., 1985).Overall our estimated C:N ratios of soil organic matter with latitude after calibration are consistent with those by Post et al. (1985).Our estimated mean C:N ratio of soil organic matter is highest for the deciduous needle leaf forest (C:N=30 g C/gN), close to the mean of the soil C:N ratio of boreal rain forest by Post et al. (1985), and much higher than those of other boreal forests by Post et al. (1985).
Assuming that the fractions of different soil P pools within the rooting zone as represented in our model are the same as those estimated for Cross and Schlesinger (1995), we estimated three parameters, v pmax for each biome (Table 1), and k plab , s pmax for each soil order (Table 2).We assumed that v pmax is a biome-dependent model parameter, because both plant roots and soil microbes can produce phosphates and growth of soil microbes depends on the supply of soluble carbon from root exudates (Treseder and Vitouske, 2001).Partitioning of soil P among labile, sorbed and strongly sorbed pools depends on soil chemical and physical properties (Barrow, 1978), and soil pedogenesis (Walker and Syers, 1976), therefore we assume that k plab and s pmax vary with soil order.62 matter by CASACNP (red) and Post et al. (1982Post et al. ( , 1985) ) (black).We calibrated 1552 using the latitudinal variation of C:N ratio of soil organic matter from Post e 1553 (Figure 4c).1554  et al. (1982, 1985) (black).We calibrated the model using the latitudinal variation of C:N ratio of soil organic matter from Post et al. (1985) (Fig. 4c).
Using nonlinear parameter optimization and the spatial distribution of IGBP biomes and soil orders, we found that the estimates of the fractions of soil organic P for different soil orders by Cross and Schelsinger (1995) do not constrain well the estimates of v pmax for all biomes.Based on the correlation of the estimates of v pmax among different biomes, we grouped the IGBP biomes into 4 groups and estimated the mean v pmax for each group (Table 1).Using the estimates of the fraction of labile, sorbed and strongly sorbed P from Cross and Schelsinger (1995) and soil order spatial distribution, we estimated p plab and s pmax for each soil order using a nonlinear optimization technique (Table 2).
Our calibrated model, with biochemical P mineralization, estimates global total soil P fractions of 34%, 9%, 10% and 47% in the soil organic matter, labile, sorbed and the strongly sorbed pools, compared to 28%, 9%, 13% and 50% from Cross and Schlesinger (1995) for the top 15 cm soil and the USDA soil order maps excluding the occluded P (Fig. 5).The difference between the two estimates is largest for Oxisol, our estimated fraction of organic P is too high and the fraction of strongly sorbed P is too low.More field measurements of biochemical P mineralization are needed, particularly for under-sampled soil orders and deeper soils (>15 cm).  Figure 5 also compares the modelled P fractions of different soil pools with or without biochemical P mineralization in the soil.Without biochemical P mineralization, the modelled fraction of P in soil organic matter accounts for over 50%, and the fractions of labile and sorbed P together are <10% for most soil orders.While the fraction of labile P in soil can vary during the growing season (Townsend et al., 2007), the fraction of P in soil organic matter is usually less than half of total P excluding occluded P for most soils (Cross and Schlesinger, 1995) except some highly weathered soil in the tropics.Consequently including biochemical P mineralization is very important for correctly representing soil P dynamics.

Steady-state pool sizes and fluxes for 1990's
Our carbon cycle is based on CASA' model (Fung et al., 2005) with some significant differences.For example, we derived the leaf phenology from the estimates of remote sensing observations (Zhang et al., 2006).Although our modeled soil C pool sizes are quite similar to those by CASA model, our modeled plant and litter pool sizes would be quite different from those by CASA'.
Estimates of carbon pool sizes at equilibrium by our model are 520 Gt C in plant biomass, 122 Gt C in litter and 2124 Gt C in soil.Overall our estimate of plant live biomass carbon (Fig. 6) shows two large peaks, one being in the tropics (15 • S to 15 • N) and the other being in the temperate and boreal region (50 • N to 65 • N).These regions account for 38% and 20% of total plant live biomass carbon.Here we compared the estimates of vegetation C pools with those by Olson et al. (1985).Because Olson et al. (1985) used different biome classification, we calculated land area weighted means of the median, minimum and maximum plant live biomass C for each 2 • latitudinal band from their spatially explicit (0.5 • by 0.5 • global) estimates.Figure 6a shows that our model vegetation biomass C agrees quite well with the mean of the median value by Olson et al. (1985) at different latitudes except two regions: the tropical region (15 • S to 15 • N) and southern temperate region (37 • S to 45 • S).
In the tropical region, where tropical forest and tropical savanna dominate, our estimated mean biomass C is 13063 g C m −2 for tropical evergreen broadleaf forest and 6220 g C m −2 for woody savanna, much higher than the mean median values of 7467 g C m −2 and 4029 g C m −2 by Olson et al. (1985) respectively.However our estimated plant live biomass carbon compares well with the estimates of 12 100 g C m −2 by Dixon et al. (1994) and 19 428 g C m −2 by Saugier et al. (2001) for tropical evergreen forest.After accounting for the difference in the area of tropical evergreen broadleaf forest used for different studies, the total plant live biomass carbon as estimated by CASACNP is 211 Gt C, similar to the estimates of 212 Gt C by Dixon et al. (1994), 244 Gt C by Ajtay et al. (1979), but much lower 65 Figure 7. Zonal mean for land grid points of fine litter production (red) and coarse woody litter production (blue) estimated by CASACNP as compared with those by Matthews (1997) (black for fine litter and grey for coarse woody debris).than the 340 Gt C by Saugier et al. (2001) for tropical evergreen broadleaf forest using the area from the IGBP vegetation map (Fig. 6b).

Latitude (degree)
In the other region between 37 • S and 45 • S, our model estimates are closer to the maximal value by Olson et al. (1985).The mean plant live biomass carbon density as estimated by CASACNP is 4605 g C m −2 , much higher than the mean median estimate of 2401 g C m −2 for this region by Olson et al. (1985).A relatively small area of land and few field measurements available may contribute to the difference between the two estimates.The region is dominated by perennial grasslands (51%) in New Zealand and Argentina (Fig. 6b) where there are few estimates of plant live biomass carbon density.
We also compared our estimates of litter productions and coarse woody debris pool sizes for different biomes with other estimates (Fig. 7); such a comparison was not previously done for the simulations by CASA or CASA'.Matthews (1997) estimated fine and woody litter production for each of 30 biome types.Using her estimates of litter production and the 1 • by 1 • biome type map of Matthews (1983), we derived the estimates of fine and woody litter productions for each 2 • latitudinal band between 60 • S to 75 • N.For CASACNP, fine litter production is calculated as the sum of litter fall from leaves and roots.
Our estimates of global fine litter production per year and the total fine litter pool size (metabolic and structural litter) are 45 Gt C year −1 and 61 Gt C, in good agreement with Matthew's (1997) estimates of 45 to 55 Gt C year −1 and 80 Gt C respectively.Our estimate of fine litter production is more variable with latitude than that of Matthews (1997), particularly in the southern hemisphere (Fig. 7).The larger fluctuation of the predicted fine litter production by CASACNP in the southern hemisphere is associated with the change in the proportion of forested land area (Fig. 6b).This regional change in biome type and the impact on fine litter production may not be estimated correctly using the empirical relationship by Matthews (1997); more field studies are needed to verify our estimates.
Estimates of woody litter production by CASACNP agree quite well with those by Matthews (1997) (Fig. 7).Our estimate of CWD flux is 6.3 Gt C year −1 and total CWD pool size is 60 Gt C globally, compared with 6.0 Gt C year −1 and 75 Gt C by Matthews (1997).Direct measurements of CWD flux are rare, as it requires successive inventories of the same plots over more than several decades, particularly in oldgrowth forests (Harmon et al., 1993).Most studies estimate the CWD production using the woody biomass and mortality rate.These estimates can be quite sensitive to infrequent disturbance, such as insect attack and extreme weather conditions.
Measurements of total CWD pool sizes are relatively straightforward and more measurements are available.Our estimates of CWD pool sizes for all forest biomes fall within the range of previous estimates.The biome mean CWD pool size we estimated is 2437 g C m −2 for evergreen needle forests, 3000 g C m −2 for deciduous needle forests, 3762 g C m −2 for the temperate and boreal mixed forests, and less than 1000 g C m −2 for tropical forests (due to rapid decomposition of woody litter in the tropics).Our estimates are comparable with the estimates compiled by Tang et al. (2003) for various forests from field measurements.The estimates they compiled vary from 1400 to 5800 g C m −2 in coniferous forests and 1380 to 2040 g C m −2 in the mixed forest in North America, and 190 to 385 g C m −2 for dry tropical forests in Venezuela, and 650 to 8500 g C m −2 in tropical rainforests in Chile, Australia and China.
Our estimate of equilibrium soil carbon of 2124 Gt C is for the entire rooting zone within which the vertical root biomass distribution is modelled using the model developed by Jackson et al. (1996), and is therefore much higher than the estimate of 1500 Gt C of Post et al. (1982) for the top 1 cm soil, but quite close to the estimate of 2300 Gt C for the top 3 m by Batjes (1996) for soil carbon.
The equilibrium nitrogen pool sizes are 6.6 Gt N for plant, 1.1 Gt N for litter and 126 Gt N for the soil organic matter, and 0.5 Gt N in the soil mineral N pool for the global terrestrial biosphere under the present climate and CASA NPP input.
Few estimates of total N in pools are available for the global terrestrial biosphere.Our estimate of total N in plant biomass is similar to the estimate of 5.6 Gt N by Xu-ri and Prentice (2008), but is much higher than the estimates of 3.1 Gt N by Gerber et al. (2010) and 3.8 Gt N by Zaehle et al. (2010).Our estimate of soil organic N is between the 100 Gt N by Post et al. (1985) for the top 1 m soil and 156 Gt N by Batjes (1996) for the top 3 m of soil globally, and quite similar to other model estimates (Zaehle et al., 2010;Gerber et al., 2010) for the rooting zones as specified in their respective models, but much higher than the estimate www.biogeosciences.net/7/2261/2010/Biogeosciences, 7, 2261-2282, 2010 of 67 Gt N by Xu-ri and Prentice (2008).Our estimate of total soil mineral N is lower than the estimate of 0.9 Gt N by Xu-ri and Prentice (2008).
There is no global estimate of total soil P for different biomes.The total amount of soil P is closely related to the property of the parent material and soil age, and the fraction of the soil P available for plant uptake is closely related to soil sorption capacity (Barrow, 1978).To estimate the amount of soil P, we used soil order to distinguish different soil mineralogy and age.Unlike the C and N cycles, most P on land is present in rocks, predominately in apatite.During pedogenesis, the phosphorus in soil parental material is mineralized into soil by weathering and uplift (Porder et al., 2007).Walker and Syers (1976) postulated that the fraction of soil P in the occluded pools unavailable to plants or soil microbes increases as soil ages.This hypothesis is supported by measurements of soil P from sites along Chronosequences in Hawaii (Crews et al., 1995) and New Zealand (Porder et al., 2007).
Estimates of global total amount of P in the terrestrial biosphere are few, and quite variable, ranging from 40 Gt P (Smil, 2000) to 200 Gt P (Jahnke, 1992).Most measurements of soil P were made on available P that only accounts for 3 to 10% of total soil P in agricultural soils, and measurements on forest soil are relatively scarce (Johnson et al., 2003).Our model estimates that the total P in soil excluding occluded P is 16.5 Gt P if biochemical P mineralization is neglected or 30.5 Gt P otherwise.Biochemical P mineralization lowers the estimate because it increases the flux from soil organic P to labile P that can be lost by leaching.Estimates of P pool sizes at equilibrium by our model are 0.40 Gt P in plant biomass, 0.04 Gt P in litter and 5.7 Gt P in soil organic matter, and 1.5 Gt P, 1.7 Gt P and 7.6 Gt P in labile, sorbed and strongly sorbed P pools in the soil.Smil (2000) pointed out that the early estimate of total soil P by Jahnke (1992) was too high, and she estimated the amount of P in soil to be 5 to 10 Gt P in organic forms and 35 to 40 Gt P in inorganic forms.Mackenzie et al. (2002) estimated that the total organic P is only about 5 Gt P globally, similar to our estimate of 5.7 Gt P in soil organic matter.Assuming that the average amount of occluded P is 35% of total soil P globally (Cross and Schlesinger, 1995), we estimate that the total amount of occluded P is 9 Gt P, and total soil P including occluded P will be 26.5 Gt P, similar to the lower estimate by Smil (2000).
Previous studies estimated that the total amount of P in terrestrial plants varies between 0.5 to 3 Gt P (Jahnke, 1992;Smil, 2000).Given the amount of N in total terrestrial plant live biomass is 6.6 Gt N, and the N:P ratio can vary from 10 to 20 g N (g P) −1 in plants (Vitousek, 1984(Vitousek, , 2004)), we conclude that the estimate of 3 Gt P in plant live biomass is too high unless we underestimate the total amount of N in plant live biomass by an order of magnitude.On the contrary we may have overestimated the amount of N in plants, as the C:N 66 side.The dotted squares represent the global terrestrial biosphere with three major 1578 compartments, plant biomass (B), litter (L) and soil (S).The units of pool size are Gt C, N 1579 or P, the mean residence time is in years and the flux units are Gt C, N or P per year.1580 Here we included occluded P in soil (see section 5.2 for further details  ratios we used are relatively low, compared with some other estimates (e.g.Vitousek, 1984Vitousek, , 2004)).Figure 8 summarizes the pool sizes and fluxes of C, N and P at steady state for the 1990's for the global terrestrial biosphere with the NPP estimates of Randerson et al. (1997) as input to our model.At steady state, the total carbon flux from plant to litter is equal to NPP, and is equal to soil respiration.The global mean NPP is 51 Gt C year −1 .The total N loss rate from soil is 0.295 Gt N year −1 , and is equal to total N input at steady state.Total plant N uptake is equal to net N mineralization, and is 1.1 Gt N year −1 , which is very close to the estimated total plant N uptake rate of 1.08 Gt N year −1 by Xu-ri and Prentice (2008).We also estimated that the annual N loss from the terrestrial biosphere is 0.06 Gt N year −1 , which is quite similar to the estimated total export of N from land to river and coastal oceans of 0.07 Gt N year −1 (Seitzinger et al., 2004).Our estimate of total N gaseous loss to atmosphere is 0.24 Gt N year −1 , and is twice as much as the global soil denitrification rate of 0.12 Gt N year −1 as estimated by Seitzinger et al. (2006).Some of the difference between the two estimates may result from N gaseous loss from nitrification and asymbiotic N fixation that is not accounted for by Seitzinger et al. (2006).
The total input of P to the terrestrial biosphere is 0.016 Gt P year −1 ; P weathering, inorganic P fertilizer addition and dust P deposition account for 12%, 84% and 4% of the total input, respectively.The rate of P loss by leaching is estimated to be 0.014 Gt P year −1 , and about 0.002 Gt P year −1 is transferred to the occluded P pools with a residence time >100 years.Using nutrient data from major rivers and coastal regions and water fluxes, Seitzinger et al. (2006) estimated the total P lost to the river and coastal ocean is 0.01 Gt P year −1 .
The mean residence at steady state can be calculated as the ratio of pool sizes and influx for C, N and P in plant, litter and soil.The total mean residence time in the terrestrial biosphere is 54 years for C, 124 years for N and 437 years for P (Fig. 8).For nutrients N and P, the exchange fluxes between plant, litter and soil within the terrestrial biosphere are much larger than the external flux into the terrestrial biosphere, therefore internal cycling of the nutrients dominates the cycling of N and P, as compared with the C cycle.The mean residence time constants of N and P in plants or litter are quite similar for N and P, but much shorter than the respective mean residence time of C, as a result of nutrient resorption by plants.

Global nutrient limitation to net primary productivity and its uncertainty
Figure 9 shows the variation of leaf N:C and P:C ratios and the nutrient limitation factor for all land points not covered by permanent snow and ice.Leaf N:C ratios of tropical forests, savannah and crop land vary between 0.04 to 0.06 g N (g C) −1 , and are significantly higher than other biomes.The N:C ratio is lowest in the deciduous needle leaf forests in the boreal region, varying between 0.02 and 0.03 g N (g C) 1 .The leaf P:C ratio varies between 0.001 and 0.003 g P (g C) −1 for unmanaged biomes, and is about 0.004 g P (g C) −1 for crop land.
Figure 9 also shows that the NPP of tropical evergreen forest and savannah and some crop land in the USA, Asia and Australia is limited by P. Most other biomes are limited by N. The deciduous needleleaf forests and high latitude shrub lands (or tundra) are most strongly limited by N.
Using Eq. ( 9) with the estimated uncertainty of leaf N:P ratio and the assumed range of C:N ratio for different biomes, we calculated the uncertainty of nutrient limitation factor (x np ) for each cell.Figure 10a shows the uncertainty of nutrient limiting factor is quite high (>0.15,shown in red) for some shrublands, grassland and woody savannah.Because of the large uncertainty of leaf N:P ratio for grassland (Table 1), we cannot distinguish N limitation from P limitation.For some shrubland and woody savannahs, the leaf N:P ratio can be quite close to 16, and therefore they are likely to be co-limited by N and P. In Fig. 10b, we show which regions are N-limited (blue-green color), P-limited (pink region) or N and P co-limited (golden color).When NPP is co-limited, the N-limiting factor is not statistically significant from the P-limiting factor (at a 95% significance level).The land type of permanent snow and ice (white) are not modeled.In order to show both N and P limitation variation spatially in the lower panel, we plotted the value of x np -1 if x n < x p , or 1-x np if x n > x p , where x n is the N limiting factor on NPP, and x p is the P limiting factor on NPP, x np = min (x n , x p ). Therefore regions with a negative value are limited by N and regions with a positive value are limited by P. A value of -0.2 corresponds to x np = x n =0.8, therefore addition of N fertilizer can increase NPP by 20%, similar for a P-limited region with a value of 0.2.
Our results agree broadly with results from a recent synthesis by LeBauer and Treseder (2008), who showed that nitrogen limitation is widespread and the relative increase in NPP in response to N fertilizer application varies from 11% for desert ecosystems to 35% in the tundra, with a global mean response of 29%.For the N-limited biomes, we estimate that N limitation reduces NPP by 10% to 40% under the present climate and (CO 2 ).LeBauer and Treseder (2008) also showed that N fertilizer addition would increase the NPP of tropical forests by about 20%, whereas our results show that nearly all tropical forests are P limited and will therefore not respond to N fertilizer addition.This discrepancy can be explained by two factors: the first one is that our model only captures the broad variations of nutrient limitation because of the relatively coarse resolution ( 2 of leaf N:P ratio and limiting nutrients are not well captured by our model simulation.For example, it has been observed that leaf N:P ratios and available soil N or P are quite variable in space and time in the tropical forests in South America (Townsend et al., 2007).The second factor is that N addition may increase biochemical P mineralization and therefore will increase NPP even when NPP is P-limited (Houlton et al., 2008).To represent this connection between N and P cycles in soil, we need to model the N cost of P uptake and N fixation explicitly.

Limitations and future studies
In this study, we developed and implemented nutrient (N and P) cycles into a well-calibrated carbon cycle model.The model was used to estimate the pool sizes and fluxes of N and P cycles, nutrient limitation and its uncertainty under the present climate condition.The estimated variation of nutrient limitation globally is also consistent with other ecological studies (Vitousek and Howarth, 1991;Le Bauer and Treseder, 2008).However, many parameters or parameterization of some processes in this model are poorly constrained, particularly those relating to nutrient cycles, because global estimates of nutrient pools and fluxes based on field measurements are very limited.To partially overcome this limitation, we assumed that all land points are at steady state under the present climate and (CO 2 ) conditions.This assumption eliminates the dependence of the modelled pools and fluxes on their initial values, and allows us to compare our model results with other estimates of pools and fluxes based on field measurements taken under present climate conditions.However the terrestrial biosphere has rarely been at steady state, particularly over the last 150 years as a result of changes in climate, atmospheric (CO 2 ), land use, and disturbance.In the following we will discuss how much the various pools have changed over the last 150 years, and what impact disturbance, such as fires and land use change, may have on the modelled nutrient limitations.
Most studies of the global carbon cycle on land assume that all pools are at steady state around 1850 (e.g.Friedlingstein et al., 2006), and integrate the model forward using prescribed inputs and simulated climate.Simulations by 11 global climate models with fully coupled carbon cycle models showed that the simulated plant biomass carbon increased by 9% and soil carbon by 4% on average as a result of changing Climate and increasing atmospheric (CO 2 ) by 2000 (Friedlingstein et al., 2006).Land use change can also affect the pool sizes and fluxes of all three cycles.It is estimated about 148 Gt C has been released into the atmosphere globally from land use change from 1850 to 2000 (Houghton, 2008) equivalent to a reduction of about 6% in the total terrestrial carbon pool we estimated at steady state.Overall the total amount of terrestrial biospheric carbon in 2000 is likely within ± 10% of the estimate for 1850 globally.
Changes in the global total size of N pools from 1850 to the 1990's will likely be small, as only 9% of the anthropogenic N input is accumulated in terrestrial pools with a residence time longer than decades (Schlesinger, 2009).However the spatial pattern of changes in carbon or nutrient pools and nutrient limitation since 1850 can be much larger.Most increases in carbon pool sizes in the terrestrial biosphere are from unmanaged forests where the only N input is deposition.As N deposition has increased since 1850, particularly in the USA and Europe, we may have underestimated the extent of nutrient limitation, particularly nitrogen limitation, for unmanaged forests.This is also consistent with the positive response of net C uptake to N deposition observed for some temperate forests (Thomas et al., 2010), since the net C uptake would be zero at steady state.For managed ecosystems, our steady state assumption will likely lead to more significant biases in the estimated nutrient limiting factors, and the biases depend on when the managed land was converted from native vegetation, and how the land was managed.For example, soil tillage may make some of the occluded P available for plant uptake while liming Can help restore soil ion balance and increase available P in the acidified soil.These have not been accounted for in our model.
Disturbance, such as fires can also have a significant effect on nutrient cycles and nutrient limitation (Certini, 2005), with the impact depending strongly on the intensity and frequency of fires (Hart et al., 2005).As Herbert et al. (2003) showed, disturbed tropical forests can become N limited, because of the disproportionately larger amount of N than P lost from timber extraction and slash burn.It is also well known that N fixers can invade after fire or during the early succession of forests in temperate regions (Vitousek et al., 2002).Therefore a steady state assumption may lead to underestimation of N limitation in those systems.In many tropical ecosystems, such as savannahs, fires can also release the phosphorus locked up in woody tissue through ash, and dispersion of ash can fertilize the vegetation regrowth after fires (Escudey et al., 2009), therefore accelerating the phosphorus cycling through different pools.Overall our steady state assumption may have resulted in overestimating phosphorus limitation in these systems; further studies are needed.
This study represents the first step in the studies of the interactions between global biogeochemical cycles and climate change.By matching the nutrient-limited NPP with the estimates of NPP by Randerson et al. (1997), we derived the estimate of nutrient limitation globally.This study has addressed the question as to what the nutrient limitation should be for the given carbon cycle at present; the question of how the carbon cycle and nutrients will interact in the future still remains unanswered.In the future, we will implement the biogeochemical model into a global land surface model that calculates the nutrient unlimited NPP as a function of a number of environmental drivers and disturbance.The combined model will then be used to study the effects of increasing (CO 2 ), land use and land use change in the past and future on pools and fluxes of all three cycles and the feedback to climate in the future.Some of the model parameters are poorly constrained, such as v pmax and values of other model parameters are arbitrarily chosen for this study.Sensitivity studies (not shown here) showed that varying v pmax by ± 20% from its mean estimate (see Table 1) only has a small influence on the estimate of nutrient limiting factors for most land points (relative change <5%).Some model parameters are assumed to vary with biome only, and will unlikely capture variations within a biome, such as the observed variation of leaf N:P ratio with species composition or seasons (Townsend et al., 2007).In the future, we will use measurements collected from field studies of ecosystem response to elevated (CO 2 ), increasing N deposition and soil warming to assess the modelled responses of different ecosystems, and therefore improve the representation of some key processes at ecosystem scale.This is important for improving our confidence in the model predictions under future climate and (CO 2 ) conditions.

Conclusions
We developed a global model of C, N and P cycles for the terrestrial biosphere.Estimates of C, N and P pool sizes and major fluxes between plant, litter and soil agree well with various independent estimates.
Including biochemical P mineralization is important for modeling the P cycle in the terrestrial ecosystem.If biochemical P mineralization is not accounted for, the model will overestimate the fraction of soil organic P and underestimate the fractions of P in the labile, sorbed and strongly sorbed pools, and the dynamics of soil P incorrectly.
Using our model for present climate conditions, we derived a spatially explicit estimate of nutrient (N and P) limitation globally that is consistent with limited evidence from field measurements.Our result shows that most tropical forest and savannahs are P-limited, and their net primary productivities are reduced by 20% due to P limitation.Most of the remaining vegetation is N-limited, and N limitation is strongest in the deciduous needle leaf forest at high northern latitudes, where N limitation reduces its NPP by about 40%.nitrogen resorption coefficient of plant pool i (=0.5 for leaf, =0.9 for wood and root) r p,i phosphorus resorption coefficient of plant pool i (=0.5 for leaf, =0.9 for wood and root) v pmax biome-specific maximal specific rate of biochemical P mineralization (d −1 ) x np nutrient limiting factor (dimensionless) x npmax a biome-dependent empirical parameter representing the ratio of nutrient un-limted NPP and nutrient-limited NPP under the present conditions (dimensionless) k n,up an empirical parameter relating plant nitrogen uptake rate to soil mineral N amount (=2 g N m −2 ) k p,up an empirical parameter relating plant P uptake rate to labile P pool size in the soil (=0.5 g P m −2 k plab an empirical parameter for describing the equilibrium between labile P and sorbed P (g P m −2 ) k ptase an empirical parameter for phosphates production (=150 g N/gP) s pmax maximum amount of sorbed P (g P m −2 ) µ i turnover rate of a plant pool i (d −1 ) µ j turnover rate of a litter pool j (d −1 ) µ k turnover rate of a soil pool k (d −1 ) µ sorb rate constant for sorbed P (d −1 ) µ ssb rate constant for strongly sorbed P (d −1 ) λ pup N cost of P uptake (=40 g N/g P for tropical biomes and 25 g N/g P for other biomes) λ ptase biome-specific N cost of phosphatase production (=15 g N/g P)

Fluxes
F c net primary productivity (g C m −2 d −1 ) F c,1990 net primary productivity for 1990's

The carbon cycle
The carbon cycle model is based on CASA' model (Fung et al., 2005) except that we combined the above ground with the below ground metabolic or structure litter pools, so that there are only nine carbon pools in our model.The equations governing the change of C pools are the same as those described by Randerson et al. (1996).They are: Where C denotes pool size in g C m −2 and µ turnover rate in d −1 , both with one subscript, i for plant, j for litter or k for soil.a c,i is the fraction of NPP (F c ) allocated to plant pool i, b j,i is the fraction of litter fall from a plant pool i allocated to litter pool j, and c k,j is the fraction of litter carbon that enters soil pool k, and d k,kk is the fraction of decomposed C from soil pool kk to soil pool k, m n is the N limitation on litter decomposition (Eq.C12), and varies from 0 to 1. Coefficient a c,j (the fraction of NPP allocated to leaf, wood or root) depends on leaf phenology.Global leaf phenology for all biomes is derived from the estimates of remote sensing observations (Zhang et al., 2006).Leaf growth is divided into four phases.Phase 1 is from leaf budburst to the start of steady leaf growth, phase 2 is from the start of steady leaf growth to the beginning of leaf senescence, phase 3 represents the period of leaf senescence and phase 4 is from the end of leaf senescence to the start of leaf bud burst.During phase 1, a c,leaf is set to 0.8, and a c,wood and a c,root are set to 0.1 for woody biomes, and 0 and 0.2 respectively for non-woody biomes.During steady leaf growth (phase 2), the allocation coefficients are constant but vary from biome to biome, taking their values from Fung et al. (2005).During phases 3 and 4, the leaf allocation is zero and its phase 2 allocation is divided between a c,wood and a c,root in proportional to their allocation coefficients.For evergreen biomes, leaf phenology remains at phase 2 throughout the year.
Leaf turnover rate will increase with cold and drought stress, and is modeled following the approach of Arora and Boer (2005).The partitioning Coefficient, b j,i , c k,j , d k,kk , µ j and µ k use the same values as in the CASA model (Randerson et al., 1996).Uptake of N by plants from soil, F n,up is calculated as F n,up = i a c,i F c n max,i − n min ,i − r n,i µ i N i (C7) N min N min + k n,up + F n,upmin where n min,i and n max,i are the minimal and maximal N:C ratios of plant pool i, k n,up is an empirical constant (=2 g N m −2 ) (Melillo et al., 1993).Here we assumed that N uptake above the minimal uptake (F n,upmin ) is proportional to the maximal amount of N required by plant growth (a c,i F c n max,i ).When the rate of plant N uptake as calculated by Eq.C7 is greater than N min /g t, (where t is the time step of model integration and equals 1 day), F n,up is set to N min / t.The minimal N uptake for a given NPP, F n,upmin , is calculated as where the first and second terms on the right-hand side of Eq. (C8) represent the minimal N required for growth and resorption, respectively.Allocation of the N uptake among different plant pools, a n,i is calculated to be proportional to the demand of N by pool i.That is = a c,i F c n max,i − n min,i − r n,i µ i N i

Nmin
Nmin+kn,up + a c,i F c n min,i − r n,i µ i N i F n,up Following Wang et al. (2007), net N mineralization rate, F n,net is calculated as the difference of gross N mineralization (F n,gr ) and N immobilization (F n,im ) rate.Gross N mineralization rate, F n,gr is calculated as The first and second terms on the right-hand side of Eq.C10 represent gross N mineralization of litter and SOM decomposition.Similarly N immobilization rate, F n,im , is calculated as m n c k,j µ j N j + kk d k,kk µ kk N kk kk =k (C11) where m n , the N-limiting factor of decomposition is calculated as follows: where F * n,net is N unlimited net mineralization, and is calculated using equations C10 and C11 with m n =1.When the N unlimited gross mineralization is less than the N unlimited immobilization, F * n,net <0 and decomposition rate is reduced by m n .A similar method was used by Comins and McMurtrie (1993) for their model.Two pathways of N loss are modeled.One is gaseous loss and the other is leaching.Gaseous N loss is proportional to net N mineralization (Firestone and Davidson, 1989) and leaching loss is proportional to the mineral N pool.That is F n,loss = f ngas max 0,F n,net + f nleach N min (C13) where f ngas is equal to 0.05 (Parton et al., 1987) and f nleach is equal to 0.5 year −1 (Hedin et al., 1995).

Appendix D Phosphorus cycle
The phosphorus cycle is based on the model of Wang et al. (2007) a p,i = 1.Where a p,i is the allocation of P uptake to different plant pools, F p,up is the plant P uptake (g P m −2 d −1 ), r p,i is the P resorption coefficient, µ sorb and µ ssb are rate constants for the sorbed and strongly sorbed P pools, respectively, both are equal to 0.0067 year −1 .s pmax and k plab are the maximum amount of sorbed P (g P m −2 ), and the constant for the adsorption (g P m −2 ), both parameters vary with soil order (Table 2).F p,net , F p,dep , F p,fert , F p,we , F p,up and F p,loss are the net biological P mineralization, dust P deposition, fertilizer P addition, P weathering rate, plant P uptake rate and P loss rates, respectively; all are in g P m −2 d −1 .F p,tase is the biochemical P mineralization rate in g P m −2 d −1 , and is calculated (see Wang et al., 2007) as F p,tase = v pmax λ pup − λ ptase λ pup − λ ptase + k ptase µ slow P slow + µ pass P pass (D11) where v pmax is the maximum specific biochemical P mineralization rate (d −1 ), λ pup and λ ptase are the N cost for P uptake and phosphatase production (g N (g P) −1 ), respectively.k ptase is an empirical constant.λ ptase =15 g N (g P) −1 and k ptase =150 g N (g P) −1 (Wang et al., 2007).Here we assumed that the contribution to biochemical P mineralization from the slow or passive SOM pool is proportional to the turnover flux of that pool.Parameters λ pup and v pmax are biome-dependent, and their values are listed in Table 1.
The dynamics of P in plant, litter and soil microbial biomass pools are similar to those for N (comparing equations D1 to D5 with C1 to C5).The treatment of the soil P pools (Eq.D8 to D10) follows that of Wang et al. (2007) except that we represent soil organic P as three pools.Different from soil N pools, the organic P in slow and passive pools can be biochemically mineralized, and the contribution to biochemical P mineralization is assumed to be proportional to their respective decay rates, As an approximation, we assumed that flux from sorbed to strongly sorbed is proportional to the size of the sorbed pool.The dynamics of the labile soil P pool is modelled using the approach developed by Lloyd et al. (2001), also used by Wang et al. (2007).
Similar to N uptake in our model, plant P uptake rate, F p,up , is calculated as F p,up = i a c,i F c p max,i − p min,i − r p,i µ i P i P lab P lab + k p,up + F p,upmin (D12) where p min,i and p max,i are the minimal and maximal P:C ratios of plant pool i, k p,up is an empirical constant (=0.5 g P m −2 ) (Wang et al., 2007).When the P uptake by plant calculated using Eq.D12 is less than P lab / t, F p,up is set to P lab / t.F p,upmin is the minimal N uptake for a given NPP and is calculated as F p,upmin = i a c,i F c p min,i − r p,i µ i P i (D13) where p min,i and p max,i are minimal and maximal P:C ratios of leaf, wood or root in g P (g C) −1 , which vary with biome type .Allocation of plant P uptake to leaf, wood and root is calculated similarly to plant N uptake.Soil P can be lost by leaching.F p,loss is calculated as F p,loss = f P P lab (D14) In this study we assumed that f P =0.04 year −1 (Hedin et al., 2003).
Fig.3.Comparison of leaf N:P (g N/g P) as estimated by CASACNP (black curve) with the empirical relationships derived from different sets of field measurements byReich and Oleksyn (2004) (dark brown),Kerkhoff et al. (2005) (orange) andHedin (2004) (yellow green).The error bars represents the one standard error of the mean leaf N:P estimate by CASACNP within each latitudinal band.The biome mean leaf N:P ratio was calibrated against the relationship of Hedin (2004).
Zonal mean of all land points of total C (a), N (b) and C:N ratio (c) of soil organic matter by CASACNP (red) and Post

Figure 5 .
Figure 5. Fraction of organic P (a), labile P (b), sorbed P (c) and strongly sorbed P (d) for each soil order excluding the occluded P for the top 15 cm soil field measurements from Cross and Schlesinger (1995) (dark red) as compared with the estimates by CASACNP with (orange) or without biochemical P mineralization (dark green).

Fig. 5 .
Fig. 5. Fraction of organic P (a), labile P (b), sorbed P (c) and strongly sorbed P (d) for each soil order excluding the occluded P for the top 15 cm soil field measurements from Cross and Schlesinger (1995) (dark red) as compared with the estimates by CASACNP with (orange) or without biochemical P mineralization (dark green).

Fig. 6 .
Fig. 6.(a) Comparison of the vegetation biomass carbon as estimated by CASACNP model (red) with those by Olson et al. (1985) (black).The grey region represents the land area weighted mean of the maximal and minimal estimates of vegetation biomass carbon, and the black curve represents the land area weighted-mean median vegetation biomass carbon as estimated by Olson et al. (1985); (b) the areas of forests, shrub land, crop land and grassland, and land ice at different latitudes.

Fig. 7 .
Fig. 7.Zonal mean for land grid points of fine litter production (red) and coarse woody litter production (blue) estimated by CASACNP as compared with those byMatthews (1997) (black for fine litter and grey for coarse woody debris).

Fig. 8 .
Fig.8.Fluxes (blue), mean residence time (red) and pool sizes (black) of the C, N and P cycles in the terrestrial biosphere at steady state under present climate conditions.The external fluxes into the terrestrial biosphere or influx are indicated on the left-hand side and the fluxes out of the terrestrial biosphere or efflux are indicated on the right-hand side.The dotted squares represent the global terrestrial biosphere with three major compartments, plant biomass (B), litter (L) and soil (S).The units of pool size are Gt C, N or P, the mean residence time is in years and the flux units are Gt C, N or P per year.Here we included occluded P in soil (see Sect. 5.2 for further details).

Fig. 9 .
Fig. 9. Spatial variation of leaf N:C ratio (g N/g C) (a), leaf P:C ratio (g P/g C) (b) and the nutrient limitation factor on NPP (c).The land type of permanent snow and ice (white) are not modeled.In order to show both N and P limitation variation spatially in the lower panel, we plotted the value of x np -1 if x n < x p , or 1-x np if x n > x p , where x n is the N limiting factor on NPP, and x p is the P limiting factor on NPP, x np = min (x n , x p ). Therefore regions with a negative value are limited by N and regions with a positive value are limited by P. A value of -0.2 corresponds to x np = x n =0.8, therefore addition of N fertilizer can increase NPP by 20%, similar for a P-limited region with a value of 0.2.

Fig. 10 .
Fig. 10.Estimates of the uncertainty (one standard deviation) of nutrient (N and P) limitation factor (x np ) (upper panel) and the regions (lower panel) where net primary productivity is N limited, P is limited or N and P co-limited.The N and P co-limited region represents the estimate of N-limitation factor (x n,leaf ) being not statistically significantly different from P limitation (x p,leaf ).
j m n µ j C j + kk d k,kk µ kk C kk − µ k C k ,kk = k;

Table 2 .
Soil order specific model parameters.Parameters k plab and s pmax are estimated in this study, P weathering rate is prescribed in this study.See 3.1 for further explanation.the uncertainty of nutrient limitation, σ xnp .Assuming that dn leaf is not correlated with dp leaf , we calculate the uncertainty of nutrient limiting factor (σ xnp ) as . There are 12 P pools.Equations governing the dynamics of P pools are:dP i dt = a p,i F p,up − µ i (1 − r p,i )P i (D1) dP str dt = (µ leaf C leaf + µ root C root )p str − m n µ str P str (D2) dP met dt = µ leaf 1 − r p,leaf P leaf + µ root 1 − r p,root P root −(µ leaf C leaf + µ root C root )p str −m n µ met P met −µ slow P slow −F p,tase µ slow P slow µ slow P slow + µ pass P pass ,kk =slow pass,j m n µ j P j + kk d pass,kk µ kk P kk −µ pass P pass − F p,tase µ pass P pass µ slow P slow +µ pass P pass ,kk = pass p,net + F p,dep + F p,fert + F p,wea + F P,tase − F p,up − F p,loss − µ sorb j c mic,j m n µ j P j + kk d mic,kk µ kk P kk (D5) −µ mic P mic ,kk = mic dP slow dt = j c slow,j m n µ j P j + kk d slow,kk µ kk P kk (D6) j c F