Implications of incorporating N cycling and N limitations on primary production in an individual-based dynamic vegetation model

Introduction Conclusions References


Introduction
A decade after the Cramer et al. (2001) seminal dynamic vegetation model (DVM) intercomparison study, the fate of the present-day terrestrial sink for atmospheric carbon dioxide (CO 2 ) in a warmer, high-CO 2 future world remains highly uncertain (Friedlingstein et al., 2006;Scholze et al., 2006;Sitch et al., 2008;Booth et al., 2012;Arora et al., 2013;Ahlström et al., 2013).Projections of change in the global net annual land-atmosphere carbon (C) flux for the coming 100 years encompass all qualitative possibilities from a markedly strengthened sink to a switch from sink to source, even if only a "business-as-usual" emissions trajectory is considered.The realised pathway varies depending on the choice of carbon cycle and climate model, the coupling between them, and other details and assumptions of a given study (Cao and Woodward, 1998;Cox et al., 2000;Cramer et al.,
One important source of this uncertainty may be the failure of many carbon cycle models to account for the constraint imposed on the production of new biomass by the plant-available pool of nitrogen (N) and its replenishment under rising [CO 2 ] (Luo et al., 2004) and climate warming (Medlyn et al., 2000).Nitrogen limits plant production in many of the world's biomes (Field, 1992), and enhanced growth in response to increasing atmospheric CO 2 or climate amelioration would require a corresponding net increase in the processes -deposition, fixation or mineralisation -governing N supply (decreased N demand or enhanced uptake of N by plants could also contribute to enabling increased growth; Finzi et al., 2007).Based on stoichiometric assumptions, Hungate et al. (2003) estimated upper and lower limits for possible future N supply and compared these with the increased demand suggested by the biomass and soil C changes projected by a suite of DVMs in the Cramer et al. (2001) study, concluding that most of the models vastly overestimated the potential increase in ecosystem carbon storage, especially the "fertilisation" effect of increased atmospheric CO 2 .In response to such criticism, ecosystem cycling of N and N constraints on plant production have been implemented in newer versions of some DVMs (Thornton et al., 2007;Sokolov et al., 2008;Xu-Ri and Prentice, 2008;Zaehle and Friend, 2010;Goll et al., 2012).Compared with C-only versions of the same models, the incorporation of N limitations generally results in lower estimates of primary production, particularly in temperature-limited ecosystems, and reduced carbon sequestration globally under future climate and [CO 2 ] projections (Thornton et al., 2007;Sokolov et al., 2008;Churkina et al., 2009;Zaehle et al., 2010;Goll et al., 2012; see Zaehle and Dalmonech, 2011 for a review on this issue).
Others have pointed to the highly generalised representation of functional diversity and spatial heterogeneity of vegetation in many current DVMs as a hindrance to their ability to accurately scale between the short-term physiological responses of individual plants to changing climate and resources, and long-term, large-scale shifts in vegetation composition and structure that underlie changes in global carbon balance (Moorcroft, 2003;2006;Purves and Pacala 2008;Fisher et al., 2010a).Most models do not keep track of size structure (demographics) of the vegetation stands they simulate and consequently fail to distinguish the contrasting resource and stress environments faced by canopy and groundlayer vegetation, which critically affect individual and population growth, succession and functional composition in closed vegetation such as forests (Purves and Pacala, 2008) and other vegetation types having a woody element, such as shrublands and savannahs.Comparing a number of DVMs widely used in carbon cycle studies to allometric data from forests around the world, Wolf et al. (2011) showed wide dis-crepancies between simulated and observed stand structure, even for sites where carbon fluxes were simulated accurately, concluding that this must inhibit the models from accurately simulating long-term changes in carbon sources and sinks.
In this paper, we present a new version of the LPJ-GUESS ecosystem model (Smith et al., 2001) that uniquely combines an individual-and patch-based representation of vegetation dynamics with plant and soil N dynamics and associated constraints on plant production and ecosystem carbon balance.We analyse the behaviour and performance of the model with respect to the main aspects of ecosystem dynamics it is designed to simulate, i.e. vegetation dynamics/biogeography, ecosystem productivity, and C and N biogeochemistry.We also revisit a number of questions previously addressed using C-only DVMs and perform new analysis to explore the implications of accounting for N cycling and its feedbacks to whole-system dynamics.These concern the enhancement of ecosystem production under elevated CO 2 and its variation among forest ecosystems globally; CO 2 and climate-driven future changes in biospheric carbon storage and the additional N supply needed to support it; and allometric scaling relationships as benchmarks of a model's ability to reproduce coupled changes in ecosystem structure and function in a consistent manner.We ask whether the conclusions drawn based on analysis with the C-only version of LPJ-GUESS remain robust or change when C-N interactions are taken into account by the model.

Ecosystem model
LPJ-GUESS (Smith et al., 2001) is a process-based dynamic vegetation-ecosystem model optimised for regional applications but also applied globally.Vegetation dynamics in the model result from growth and competition for light, space and soil resources among woody plant individuals and a herbaceous understorey in each of a number of replicate patches for each simulated site or grid cell.The suite of simulated patches is intended to represent the distribution within a landscape representative of the grid cell as a whole of vegetation stands with different histories of disturbance and stand development (succession).Individuals for woody plant functional types (PFTs; trees and shrubs) are identical within a cohort (age/size class) and patch.Photosynthesis, respiration, stomatal conductance and phenology (leaves and fine roots) are simulated on a daily time step.The net primary production (NPP) accrued at the end of each simulation year is allocated to leaves, fine roots and, for woody PFTs, sapwood, following a set of prescribed allometric relationships for each PFT, resulting in height, diameter and biomass growth.Population dynamics (establishment and mortality) are represented as stochastic processes, influenced by current resource status, demography and the life-history characteristics of each PFT (Hickler et al., 2004;Wramneby et al., 2008).
An overview of LPJ-GUESS is provided in Appendix B. The version of LPJ-GUESS used in the present study incorporates N cycling in vegetation and soil, and N limitations on plant production.It is further developed from the C-only standard version of the model, version 2.1, which is fully described in Ahlström et al. (2012) and references therein.An overview of the N-cycling scheme follows; further details are provided in Appendix C.
A pool of mineral N -available to plants, via root uptake, and to soil microbes -is provided by atmospheric deposition, biological N fixation and gross N mineralisation of soil organic matter (SOM).Annual N deposition was prescribed externally to the model (see below).Biological N fixation is estimated as a dynamic function of evapotranspiration (simulated by the model), following a parameterisation proposed by Cleveland et al. (1999) that builds on an empirical correlation between N fixation and evapotranspiration at the global scale.
C and N dynamics of soils are simulated conjointly by an SOM scheme adopted from the CENTURY model (Parton et al., 1993), with modifications by Comins and McMurtrie (1993) and Kirschbaum and Paul (2002), and recent updates by Parton et al. (2010).Decomposition of 11 SOM compartments differing in C : N stoichiometry and resistance to decay results in respiration (release of CO 2 ) and transfer of C and N between pools, satisfying mass balance.Carbon entering the receiver pool drives N mineralisation or immobilisation.Gross N mineralisation is effected as an increase in the soil mineral N pool when N transferred from a donor pool exceeds the corresponding increase dictated by the C : N ratio of a receiver pool (N "supply" exceeds "demand"); in the opposite case (N demand exceeds supply), N immobilisation occurs, reducing mineral N. Net N mineralisation, describing the net release of N from decaying organic matter, is the balance between these two terms.Decay rates are sensitive to soil temperature and moisture, and C : N ratios for certain SOM compartments change within limits depending on available N, emulating dominance shifts between decomposer functional groups as soil nutrient status changes (Rousk and Bååth, 2007).
Plant N uptake, computed on a daily time step, is the smaller component of the supply provided by the soil mineral N pool, and that part of daily N demand for allocation to the leaf photosynthetic apparatus and to biomass growth that cannot be satisfied by retranslocation from storage and from leaves, fine roots and sapwood undergoing phenological turnover.N demand is driven by optimal leaf N content, computed following Haxeltine and Prentice (1996a) as a linear function of the carboxylation capacity estimated to maximise the net photosynthesis of the plant canopy given current plant size, phenology, architecture and microenvironment.N demand for allocation to growth in other tissues follow leaf N, conserving relative differences between leaves, fine roots and sapwood in N concentration of new growth (Friend et al., 1997;White et al., 2000;Zaehle and Friend 2010).Plants retain half of the N content of shed roots and leaves and sapwood on conversion to heartwood for retranslocation to remaining tissues.The assumption of a 50 % resorption efficiency for N is consistent with reviews by Aerts (1996) and Vergutz et al. (2012).Excess nitrogen is retained in the nitrogen store, which serves to buffer the effects of seasonal and to some extent interannual variations in the balance between N demand and supply.
N limitation occurs if the N demand of vegetation in a patch, net of retranslocation, cannot be met by the supply of mineral N in the soil.Under N limitation, plant individuals take up soil N in proportion to their fine root surface area, and photosynthesis and allocation for the current year are reduced.N limitation may result in an increased relative allocation to fine roots, promoting more efficient N uptake the following year.
To diagnose and illustrate the consequences of incorporating an active N cycle in LPJ-GUESS, for a number of the model experiments described in the next section, parallel simulations were performed with N constraints on primary production enabled and disabled.These are denoted C-N and C-only, respectively, throughout the paper.All simulations were performed with the enhanced version of the model, including the CENTURY-based SOM scheme, and other modifications, as described above and in Appendix C. N constraints were disabled in the C-only simulations by artificially satisfying plant demand for N, calculated on the basis of optimal carboxylation capacity, irrespective of the available supply of mineral N in the soil; see Appendix C for further details.
With the exception of the simulation used to produce the N-limitation map (further detailed below), C-only simulations were performed with a different value of a global parameter, α a , that linearly scales the quantum efficiency of photosynthesis, equally for all vegetation and at all points in time and space (see Appendix B for a further explanation of this parameter).The parameter was calibrated to attain comparable overall levels of the main ecosystem stocks and fluxes of carbon at global scale under 1961-1990 forcing, with and without N constraints enabled.

Global hindcast experiment
To relate the performance of the model to other studies and observational benchmarks on recent historical carbon balance at global scale, and to analyse the implications of accounting for C-N interactions on the model predictions, 20th-century hindcast simulations of potential natural vegetation, C and N dynamics were performed across a 0.5 • × 0.5 • (latitude-longitude) grid covering the ice-free global land surface.Separate simulations were performed with and without N constraints enabled.Climate forcing (mean monthly temperature, precipitation and cloud fraction) was provided by the Climate Research Unit (CRU) TS 3.0 observed climate data set (Mitchell and Jones, 2005).Atmospheric CO 2 concentrations, varying annually, were prescribed from observations (McGuire et al., 2001).The simulation was initialised with a 500-year "spin-up" to build up vegetation and soil C and N pools in an approximate steady state with early-20th-century climate.CRU climate data for 1901-1930, cycled repeatedly and with any interannual temperature trend removed, were used to force the model during the spin-up.A CO 2 concentration of 296 ppmv (1901 value) was assumed throughout the spin-up.
Monthly N-deposition rates, varying decadally and across the global grid, were taken from the ACCMIP database of Lamarque et al. (2013).The values were interpolated bilinearly from their original 1.9 • × 2.5°grid to the 0.5 • × 0.5 • grid of the climate data.The 1850 deposition value for each grid cell was used up to and including the simulation year corresponding to 1850 in the spin-up.

Forest ecosystem measurement sites
To assess in more detail interactions between stand structure and ecosystem function at the local scale at which such interactions occur, site-scale simulations were performed for European forest ecosystems for which measurements of forest structure, production and N uptake are available from the Carbon and Nitrogen Cycling in Forest Ecosystems (CANIF) data set (Schulze, 2000).A subset of 10 sites, 5 dominated by evergreen conifers and 5 by winter deciduous broadleaved species, for which site climate data were available, were chosen for analysis.For each site, meteorological forcing data (monthly temperature, precipitation and incoming shortwave radiation) were constructed by regressing monthly values from the CRU data set on corresponding values from the available climate time series for the specific site, for the years for which the CRU and site data coincide (2-8 years depending on site).The resulting regression coefficients were used to construct a "site-adjusted" version of the CRU data to force the model through a 500-year spin-up and historical simulation from 1901 up to and including the period of available CANIF data.Prior to 1901, the 1901-1930 data, with any temperature trend removed, were cycled repeatedly.Site-specific N-deposition data were constructed in a corresponding manner to the climate data, based on the Lamarque et al. (2013) data for the nearest grid point and N-deposition data for each CANIF site.N deposition and atmospheric [CO 2 ] during the spin-up followed the same approach as in the global hindcast simulation.The PFT mixture simulated natively by the model (i.e.potential natural vegetation in equilibrium with the forcing climate) was retained during the spin-up.For the historical period of the simulation, actual CO 2 concentrations and N-deposition values from the nearest grid cell were applied.Forest management history (where available) was emulated by transferring all vegetation to litter in the historically recorded year of planting and establishing a new stand with the PFT composition and density prescribed from data, where available.For sites for which planting density was not known, the value was calibrated to attain the observed density of the measurement years.
As the CANIF data set is restricted to a limited range of European forest ecosystems, NPP predictions from the model were also compared to productivity data from 132 sites across the wooded biomes of the world, compiled by Luyssaert et al. (2007).Simulation data from the nearest 0.5 • × 0.5 • grid cell for 1961-1990 from the global hindcast experiment (above) were extracted for comparison to the Luyssaert data.

Forest succession experiment
To explore and demonstrate the ability of LPJ-GUESS to reproduce broad features of the structural and compositional evolution of forest stands during secondary succession (i.e.recovery following a disturbance event), observed data on species composition and size structure for 12 virgin boreal forest stands in northern Sweden compiled by Linder et al. (1997) were aggregated into a chronosequence, i.e. rank ordered in terms of stand age.For comparison to these data, a model simulation was performed for the 9 nearest 0.5 • × 0.5 • grid points from the CRU climate database surrounding each observed stand.Following a standard spin-up (see above) all vegetation was removed, emulating a disturbance, and the simulation continued for 270 years forced by climate data for 1961-1990, cycled repeatedly.A 1901 CO 2 concentration of 296 ppmv was used to force the model initially, succeeded by the CO 2 concentration time series for  for the last 90 years of the simulation.N-deposition data from the same grid points as the climate data were employed (see Sect. 2.2.1), forcing the model with the 1850 value initially, thereafter the observed time series for 1860-1990.Wolf et al. (2011) examined allometric and size structure relationships from three databases encompassing biomass and stand structure data for forest ecosystems globally, comparing observed relationships to results from eight DVMs.Allometric scaling metrics, such as the relationship between mean tree mass and population density (expected to be constrained by "self-thinning" under crowding in natural stands; Enquist et al., 1998) or the relationship between stem and foliage mass (influenced for example by the trade-off between height growth and light interception during the development of a stand), reflect the relative evolution of structural and functional parameters during the stand development following clearcutting-planting or a natural disturbance.The ability of models to reproduce scaling metrics observed in real forest stands may therefore be seen as a relevant benchmark of performance in simulating conjoint changes in structure (e.g.biomass) and function (e.g.NPP) for wooded ecosystems.With one exception, the DVMs considered by Wolf et al. (2011) all employ large-area parameterisations of biomass accumulation that, in contrast to LPJ-GUESS, do not distinguish age/size classes of trees and thereby cannot discriminate population from individual growth.To examine allometric scaling in the vegetation structure predicted by LPJ-GUESS we extracted data for the calendar year 2001 from the results of the global hindcast experiment (above) for those grid cells that encompassed one or more stands from the Cannell (1982) or Usoltsev (2001) databases considered by Wolf et al. (2011).Following Wolf et al. (2011), straightline relationships were fitted by simple linear regression to log-log plots of mean tree biomass (M, kg) versus population density (individuals ha −1 ) and foliage versus stem biomass (kg), for this sample of grid cells.The resulting curves were compared with the observed stand data and the corresponding curves for the DVMs considered by Wolf et al. (2011).

Forest FACE experiment
In a study spanning forest ecosystems globally, Hickler et al. (2008) simulated the potential effect on NPP of raising CO 2 concentrations by 200 ppm in a manner emulating CO 2 treatment in free-air CO 2 enrichment (FACE) experiments.Elevated CO 2 concentrations were simulated to enhance NPP by 15-35 % depending on the biome.However, the authors conjectured that these findings might be overestimates, because nutrient limitations, not represented in the C-only version of LPJ-GUESS they used for their simulations, might limit the NPP enhancement under elevated CO 2 , as widely discussed in the literature (e.g.Luo et al., 2004) and demonstrated empirically in some elevated-CO 2 experiments (Hyvönen et al., 2007;Norby et al., 2010).
We repeated the same model experiment performed by Hickler et al. (2008) with N cycling enabled.The global model experiment was carried out using the gridded CRU TS 3.0 global climate data set (Mitchell and Jones, 2005).Two simulations were performed: one with actual historical CO 2 concentrations from 1901 to 2002(McGuire et al., 2001)), and one with historical CO 2 concentrations increased to 550 ppmv during 1996-2002.The chosen time period corresponds approximately to that in which the FACE measurements considered by Hickler et al. (2008) were obtained.
The model was run with potential natural vegetation and for all grid cells classified as forest in the Haxeltine and Prentice (1996b) map of global potential natural vegetation, with the exception of tropical deciduous forests, which have a savannah-like structure.Following a 500-year spin-up to establish the "steady-state" vegetation, the model was driven by the observed climate from 1901 to 2002.Hungate et al. (2003) estimated upper and lower limits for possible global N supply to terrestrial ecosystems over the 21st century and compared these with the increased N demand suggested by the biomass and soil C changes projected by a suite of six DVMs in a future-climate study by Cramer et al. (2001).They concluded that a majority of the models vastly overestimated the potential increase in ecosystem carbon storage, especially the fertilisation effect of increased atmospheric CO 2 , when accounting for limitations in the amount of N available for the production of additional biomass.To evaluate and compare the performance of LPJ-GUESS with and without N constraints under future-climate and [CO 2 ] forcing, and to compare the results in terms of global N demand with other models and the limits proposed by Hungate et al. (2003), we extended the global hindcast experiment, described above, to 2100 forced by climate and [CO 2 ] data from the RCP 8.5 climate scenario with the MPI-ESM-LR general circulation model (GCM; Giorgetta et al., 2013) and N-deposition data from Lamarque et al. (2011).

Global C and N stocks and fluxes
The main stocks and fluxes of the biospheric C and N cycles as simulated in the global hindcast experiment are shown in Table 1, where the results may be compared between C-only and C-N simulations, and with literature estimates based on observations and other DVM-based studies.In general, simulated values are within the ranges suggested by other studies.Exceptions include biomass burning, for which the global mean annual flux of around 5 Pg C from the C-N simulation exceeds the upper-range literature estimates by some 30 % (Ito and Penner, 2004;van der Werf et al., 2010).Fire resistance, a PFT-specific parameter, may be arguably set too low in LPJ-GUESS, resulting in a higher combustion efficiency per unit burnt area than observed in reality, e.g. for tropical evergreen trees for which 10-30 % resistance to burning (70-90 % combustion) is assumed by the model (Sitch et al., 2003), whereas in reality only about half of standing biomass may be destroyed in an average fire (Pinard et al., 1999).Simulated global N fixation is much smaller than the values suggested by observation-based studies (Table 1).N fixation in LPJ-GUESS is computed based on the "conservative" dependency on evapotranspiration proposed by Cleveland et al. (1999) (see Appendix C).The latter yields a global N fixation estimate of 100 Tg N yr −1 based on the evapotranspiration data used in the Cleveland et al. (1999) study, but the LPJ-GUESS estimate is additionally affected by a downward discrepancy of approximately 50 % in simulated global evapotranspiration (352 mm yr −1 ) compared with the evapotranspiration values, totalling 717 mm yr −1 globally, on which  In general, accounting for C-N interactions appears to have a minor influence on simulated stocks and fluxes of carbon under the present-day climate at global scale in LPJ-GUESS; however, it should be borne in mind that the model has been calibrated (by adjustment of the quantum efficiency scalar, α a ; see Appendix B), independently for the C-only and C-N simulations, to the literature estimates shown in Table 1.
The global distribution of ecosystem NPP simulated by LPJ-GUESS C-N in the global hindcast experiment is shown in Fig. 1a.The geographic pattern is comparable to that suggested by other model-based studies (Schloss et al., 1999;Running et al., 2004) and inferred from observations (e.g. for gross primary production (GPP) by Jung et al., 2011), with the highest productivity simulated to occur in the wet tropics, intermediate levels in mesic parts of the middle latitudes, declining to negligible levels in the deserts and high Arctic.
In an intercomparison study involving several DVMs, Piao et al. (2013) found that the C-only standard version of LPJ-GUESS exhibited a negative bias in GPP in the wet tropics compensated in part by a positive bias in high northern latitudes, in comparison to global GPP inferred based on ecosystem flux data by Jung et al. (2011).Mean global GPP (113 ± 3 Pg C yr −1 ) was comparable to the Jung et al. (2011) estimate of 118 ± 1 Pg C yr −1 .In the C-N simulation, N lim-itation of production associated with low mineralisation and fixation rates results in lower productivity in cold-climate areas (Fig. 1b) and a steeper gradient from high to low latitudes compared with the C-only simulation.Due to recalibration of α a , global GPP and NPP remain relatively similar to the Conly simulation (Table 1).
Apart from cold-climate areas such as the Boreal Zone and Arctic, N limitation is most strongly expressed in waterlimited ecosystems, for example desert and shrublands of western North America and the steppe belt of Central Asia (Fig. 1b).Soil moisture status affects decomposition rate and net N mineralisation in the model (Appendix C).However, the primary, decomposition-driven effect is amplified by a vegetation-dynamics-mediated secondary mechanism in which a differential degree of N limitation causes a shift from the relatively productive but N-demanding C 3 grasses towards an increased representation of woody PFTs.In the model, woody PFTs have a lower productivity under the same driving conditions due to the costs of maintaining support and transport structures (sapwood) that reduces the amount of carbon available for production of foliage.C 4 grasses, which are constrained to grow in areas with a coldest-month mean temperature of at least 15 • C, are not affected by N limitation (due to much lower carboxylation capacity, and accordingly C : N ratio, compared with C 3 grasses; Collatz et al., 1992), resulting in a stark contrast in N limitation across the 15 • C coldest-month isotherm, most obviously in northern Australia.

Potential natural vegetation
Vegetation patterns simulated by LPJ-GUESS in the C-N simulation (Fig. 2b) are comparable to results obtained by Hickler et al. (2006) using the C-only standard version of the model run in the more generalised population mode (similar to LPJ-DGVM; Sitch et al., 2003).The main differences concern a lower herbaceous component at the expense of woody vegetation in water-limited ecosystems in the present study, resulting in an increased coverage of woodlands and shrublands and less savannah and grassland.A further difference concerns the extent of the Larix-dominated boreal deciduous forest in Siberia (see below).
Comparison with the Haxeltine and Prentice (1996b) map of potential natural vegetation (Fig. 2a) reveals broad agreement in the geographic locations and ranges of major forest types versus water-limited vegetation types such as savannahs, shrublands and grasslands, as well as tundra.Agreement at the level of individual vegetation classes is poorer, but is influenced by thresholds in the classification scheme for transforming simulated PFT abundances to vegetation classes for mapping (Appendix A, Table A1).The boreal forest tree line in the Northern Hemisphere (NH) is correctly placed by the model, with the exception of eastern Siberia, where the model predicts a more southerly transition zone  A1).
to Arctic tundra, compared with the vegetation map.Tundra is also simulated for parts of eastern Alaska and northwestern Canada classified as boreal evergreen forest/woodland in the vegetation map.The C-only simulation (Fig. 2c) reproduces the observed distribution of forest versus tundra for this area, revealing that N limitation is acting too strongly in the C-N simulation, suppressing the tree line by reducing tree productivity below the limits for continuous forest cover.A more limited extent of the Siberian Larix-dominated boreal deciduous forest belt was predicted by LPJ-GUESS C-N compared both with the vegetation map and the Hickler et al. (2006) results.This cannot be explained by the coldseason temperature limits prescribed for boreal deciduous trees -which are in common with the Hickler et al. (2006) study -but reflects the influence of N limitation, reducing tree productivity and lowering leaf area index (LAI) below the threshold of 0.5 for boreal deciduous forest/woodland in the classification scheme.This is confirmed by the results from the C-only simulation (Fig. 2c), boreal deciduous forest becoming much more extensive when N limitation is switched off.
The Boreal-temperate forest transition is generally well simulated, although an arguably low threshold of 0.5 for the combined LAI of boreal needleleaved and summergreen trees (Appendix A; Table A1) may exaggerate the extent of boreal evergreen forest/woodland in Asia and western North America.N limitation reduces the southerly extent of the Boreal Forest into dry-climate parts of North America and Central Asia.Temperate forest of both hemispheres is correctly placed, although in the NH a comparatively low representation of deciduous broadleaved trees results in an overrepresentation of the classes temperate/boreal mixed forest and temperate broadleaved evergreen forest, when compared with the vegetation map.Water-limited (xeric) vegetation belts of western North America, Central Asia and the Mediterranean are predicted by the model; however, woody PFTs are overrepresented in comparison to grasses, with the result that some steppe (Central Asia) and prairie (North America) areas are misclassified as xeric woodland/shrubland instead of grassland.This apparent bias favouring woody PFTs in drier climates may be explained in part by the arbitrary LAI thresholds that determine whether a certain mixture of trees and grasses is classified as woodland/shrubland, grassland or savannah in the classification scheme (Table A1); presence of both trees and grasses is simulated in most areas mapped as xeric woodland/shrubland in the C-N simulation.However, a general shift towards a greater representation of woody vegetation, compared with the Hickler et al. (2006) study, is apparent both in the C-only and C-N simulation results of the present study.A likely explanation is an update of a parameterisation in the model expressing specific leaf area (SLA, i.e. leaf area per unit leaf mass) as a function of leaf longevity (see Appendix C).As a result of this update, SLA is reduced for all PFTs, but PFTs having shorter-lived leaves, including deciduous grasses, experience a proportionately greater reduction in PAR uptake capacity, and therefore NPP.The result is a competitive shift favouring woody PFTs at the expense of grasses, which is most strongly expressed in dry-climate areas such as savannah and steppe areas.
Equatorial rainforest of Central America and the Amazon, the Congo Basin and southeast Asia is correctly placed by the model.A high representation of tropical evergreen trees results in some areas classified as tropical seasonal forest in the vegetation map being classified as tropical rainforest by the model.For the same reason, moist savannah on the vegetation map tends to be classified as forest by the model.

Forest ecosystem measurement sites
Model performance in comparison to observations from the CANIF forest measurement sites and biome averages from the Luyssaert et al. (2007) database is shown in Fig. 3.In general, the model appears successful in simulating the canopy height (Fig. 3b) for a given stand density (Fig. 3a), bearing in mind that density was calibrated to the observed data for each site (see Methods), while height is a free diagnostic variable in the model.In the case of NPP (Fig. 3c), the average simulated for all sites by the model (0.52 kg C m −2 yr −1 ) may be compared to the average for the observations (0.43 kg C m −2 yr −1 ), but there is little agreement between the model and observations among sites.Sim-ilarly poor agreement between simulated and observed NPP for the CANIF sites was noted for the O-CN model by Zaehle and Friend (2010).When compared with biome-specific averages from the larger Luyssart database, there is moderate agreement: the model, like the Luyssaert estimates, generally exhibiting higher productivities for closed forest ecosystems of the mid-latitudes and tropics, and lower values for strongly water-and temperature-limited ecosystems such as savannahs, grasslands and some types of boreal forest.For nitrogen-use efficiency (NUE), i.e.NPP per unit vegetation N uptake (Fig. 3d), the model results span a much wider range among sites compared to the CANIF data.Vegetation and litter C : N ratios in LPJ-GUESS are assumed to vary proportionately (on allocation) with leaf C : N, the latter in turn depending on the carboxylation capacity that maximises current net C assimilation subject to the amount of N needed for production of the carboxylation apparatus (largely the enzyme rubisco) (Appendix C).This means that tissue C : N stoichiometry is allowed to vary within wide bounds depending on dynamic variations in plant demand versus supply of N. It is possible that tissue C : N ratios in the model are thus insufficiently constrained, exaggerating the structural and biochemical plasticity of real plants, and overestimating variation in NUE.Results from an intercomparison of 11 ecosystem models, including LPJ-GUESS, with N cycling effects on production enhancement inferred from the results of forest FACE experiments suggest that LPJ-GUESS and other models simulating flexible stoichiometry generally overestimate the stoichiometric plasticity of trees (Zaehle et al., 2014).Smith et al. (2001) and Hickler et al. (2004) demonstrated the ability of the C-only standard version of LPJ-GUESS to reproduce observed patterns of tree species replacement during secondary succession in temperate and boreal forest.A replacement series from fast-growing, shade-intolerant pioneer species to slow-growing, shade-tolerant and longlived species is predicted by classical successional theory and emerges dynamically in models, like LPJ-GUESS, implementing the gap-phase concept of forest dynamics (for a review see Bugmann 2001).As expected, the model exhibits similar qualitative behaviour in a C-N simulation of the recovery of boreal forest vegetation in Sweden after a complete disturbance (Fig. 4b).Pioneer species, for this region mainly birch (Betula spp.) and Scots pine (Pinus sylvestris) (Linder et al., 1997), dominate the initial tree layers but suffer reduced regeneration and/or mortality under shading.They are eventually replaced by shade-tolerant species, here Norway spruce (Picea abies), that are able to regenerate in small treefall gaps and survive, continuing to grow slowly, even in the deep shade of the closed forest canopy.In the absence of stand-replacing disturbance, pioneers are no longer regenerating and consequently absent from the smallest size   category (< 5 cm stem diameter) after about one century, persisting for longer in the canopy (larger tree size classes, > 15 cm) where shading is less pronounced.A chronosequence constructed by ordering 12 forest sites in order of age since the last disturbance suggests qualitatively similar successional dynamics and size structure compared with the simulation results (Fig. 4a).It should be noted that each point (stand age) in the chronosequence represents a unique site, so that the depicted variability reflects differences in population structure both in space (local edaphic and climatic environment, stand composition, small-scale disturbances) and time (succession).By contrast, the model results are averaged among grid cells, depicting only variability due to stand development (and a small component due to stochastic processes in the model).

Forest stand allometry experiment
In Fig. 5, allometric scaling relationships fitted to simulated data for global forests in the forest stand allometry experiment are compared to observations (all panels) and to the corresponding curves relating tree foliage mass to stem mass (Fig. 5c and d) for the eight models considered by Wolf et al. (2011).In general, the slope coefficient (β) of the depicted allometric relationship shows similar or better agreement to the observations when N limitation is taken into account (Table 2).For broadleaved forest (Fig. 5a), simulated stand density declines in a similar proportion to mean tree size (slope) as observed, but the simulated stands are generally sparser (lower intercept) than the observed forests.For needleleaves, mean density is captured more accurately.The relationship between foliage mass and stem mass may be seen as a relevant indicator of the relative evolution of function (tree light interception, NPP) and structure (tree biomass) as stands develop.For this relationship, LPJ-GUESS C-N exhibits a goodness of fit comparable to the best of the eight models considered by Wolf et al. (2011) (Fig. 5c and d).An improvement relative to the C-only simulation is apparent for broadleaves in Fig. 5c, where N cycle dynamics result in a flatter slope of the mean relationship between foliar and stem biomass (Table 2).For needleleaves N cycling increases the slope, but the fit to the observed data is closer.The dynamics of stand structure in LPJ-GUESS are an emergent outcome of combined effects of individual allometric growth, population dynamics, community-level (among PFT) interactions  and disturbance, making any improvements in accuracy difficult to trace to any particular process or interaction in the model.

Forest FACE experiment
FACE treatment applied synthetically to wooded ecosystems across the globe resulted in a qualitatively similar global pattern in a C-N simulation compared with Hickler's et al. (2008) analysis based on the C-only standard version of LPJ-GUESS (Fig. 6).CO 2 fertilisation results in an increased productivity in all regions, but the percentage enhancement of NPP is lowest in high-latitude ecosystems, increasing along a temperature gradient to the tropics.The main mechanism underlying this pattern is the greater relative enhancement of leaf-level net C assimilation at higher temperatures expected as a result of the suppression of photorespiration at higher temperatures in C 3 vegetation (Long 1991).Such a response is encapsulated in the Farquhar-based photosynthesis model adopted by LPJ-GUESS (Hickler et al., 2008).Quantitatively, however, the gradient in NPP en A further qualitative difference is the appearance of NPP enhancement in the range 40-50 % in warm temperate to subtropical parts of southeastern North America and Asia.The C-only model simulates a more moderate enhancement of 20-40 % for these areas (Fig. 6c).Such contrasting behaviour between the models may be traced to a synergistic effect on photosynthesis and autotrophic respiration, resulting in an increase in C-use efficiency, i.e.NPP as a proportion of GPP, for the temperate evergreen trees simulated in these areas in the C-N simulation.The enhanced photosynthesis resulting from the step increase in CO 2 concentration (emulating FACE treatment) raises the C content of new biomass relative to N. Since respiration scales with tissue N content, this results in a decrease in respiration that augments the enhancement of NPP, relative to a rise in GPP The same mechanism can be expected to operate for all PFTs and climate zones, but is suppressed in areas experiencing greater N limitation, such as cool temperate and boreal forests, and in areas with higher initial LAI.The latter explains why NPP enhancement is generally lower in high-LAI tropical ecosystems, e.g. in the Amazon and Congo Basin (Fig. 6c).
Whereas the majority of forest FACE experiments have been implemented in temperate forests, there is consensus based on theory and results from whole-tree chamber manipulations and one tree line FACE experiment that coldclimate ecosystems may experience much lower rates of NPP enhancement under elevated CO 2 compared with the 21-25 % increase typically seen in temperate forest FACE experiments, at least during the first years of treatment (Norby et al., 2005;Kostiainen et al., 2009;Dawes et al., 2011).Nitrogen limitation, allocation shifts from above-to below-ground parts of the plant as well as sink limitation for the additional assimilates resulting from a CO 2 -driven increase in photosynthesis are offered as explanations (Hyvönen et al., 2007).
It should be noted that the above-mentioned results concern the modelled influence of elevated CO 2 concentrations on NPP.We do not analyse the fate of the increased C uptake within the ecosystem.Results from FACE experiments and other CO 2 manipulations reveal that increased photosynthesis or NPP will not necessarily translate into larger biomass stocks or ecosystem carbon storage, for example if the extra C is mainly allocated to pools with fast turnover rates, such as fine roots or root exudates (Norby et al., 2002;Körner et al., 2005;Finzi et al., 2007).

Future-climate experiment
Under future-climate and [CO 2 ] forcing, both the C-N and C-only simulations exhibited an increase in ecosystem C storage globally of the order of 10 % (Fig. 7; Table 1).The additional N required to match this increase in C sinks falls within the availability estimates of Hungate et al. (2003) for both simulations (Fig. 7).Under [CO 2 ] forcing only, the increase in C storage is larger by a factor of 3 in the C-only simulation, but only 20 % in the C-N simulation, pointing to a marked suppression of net CO 2 fertilisation globally when N dynamics are introduced.A negative climate influence on plant and soil C was likewise simulated by all models in the Hungate et al. (2003)    The additional N requirement simulated in the C-only simulation under [CO 2 ] forcing alone only slightly exceeds the high N supply limit of Hungate et al. (2003) and is lower compared to any of the six models considered in that study.All of the latter models are "first-generation" DVMs that employ large-area parameterisations of vegetation dynamics and whose ability to replicate the temporal evolution of vegetation structure and biomass accumulation under rapidly changing driving conditions may be questioned (Purves and Pacala, 2008;Fisher et al., 2010a).In the LPJ model, for example, the ability for vegetation C pools to increase to absorb a CO 2 -driven increase in primary production is limited by a strict geometric constraint imposed by the prescribed allometry of "average individual" plants (Sitch et al., 2003).As a result, increased productivity may too strongly translate into increased biomass turnover (through phenology and mortality) and soil C storage, requiring a steady supply of N for the continuous replacement of plant biomass.
In contrast to many first-generation DVMs, LPJ-GUESS mechanistically accounts for the expected demographymediated lag in the response of biomass accumulation to whole-ecosystem NPP in areas newly rendered climatically suitable for the growth of trees, such as tundra distal of the Boreal Forest tree line under a warming Arctic climate  (Miller and Smith, 2012).Demographic lags may be expected, and are simulated by the model, even in areas already occupied by trees but where increased productivity enables stand densification to occur; such effects will be simulated under CO 2 fertilisation alone, in the absence of climate change.In addition, Piao et al. (2013) show that LPJ-GUESS simulates a lower baseline (present-day) GPP globally than current versions of all but one of the models (SDGVM) analysed by Hungate et al. (2003).This may reflect a low bias for GPP in the wet tropics in the C-only version of LPJ-GUESS (Piao et al., 2013), although simulated GPP in that study was consistent globally with an independent data-driven estimate (Jung et al., 2011;Piao et al., 2013).
In the C-N simulation, the additional N requirement falls between the low and high supply limits of Hungate et al. (2003).It may be argued that these limits are too low as they do not account for the influence on N supply of increased mineralisation in warming soils, nor for a possible increase in plant N uptake capacity at elevated CO 2 (Finzi et al., 2007).However, accepting the (conservative) logic of the Hungate et al. (2003) analysis, the estimates of C and N sequestration under future-climate and [CO 2 ] forcing suggested by the C-N simulation appear to be realistic and mutually consistent.

Discussion
Accounting for N cycle dynamics had two main overall effects on the performance and behaviour of our model.Firstly, negative effects of low temperatures and soil moisture deficits on SOM dynamics and N mineralisation resulted in a soil N-mediated relative reduction in plant productivity most strongly expressed in cold-climate ecosystems of the Boreal Zone and Arctic, and dry-climate regions of middle latitudes.Forest ecosystems occupying well-watered parts of the temperate zone and tropics were less affected by N limitation.As a result, NPP appears to exhibit stronger large-scale geographic gradients in the C-N version of LPJ-GUESS, compared with the standard C-only version of the model.This is likely to represent an improvement in performance when compared with observation-based estimates of largescale productivity (Piao et al., 2013).Secondly, N limitation resulted in shifts in vegetation structure, generally favouring woody PFTs and C 4 grasses and disfavouring C 3 grasses.This effect was most pronounced in temperate dry-climate regions such as western North America and Central Asia and would tend to amplify the primary (physiology-driven) N limitation in the model, since woody PFTs, having higher C costs due to the production of stems, exhibit lower growth efficiency (NPP per unit leaf area) than grasses when growing under the same conditions.In areas where N limitation results in a marked reduction in whole-ecosystem NPP, however, an opposing shift, favouring grasses at the expense of trees, arises from the release of grasses from suppression due to shading and preemption of space by trees.In comparison to a map of potential natural vegetation of the world, accounting for N cycle dynamics thus appeared to result in improved model agreement in some areas (e.g.western North America, Central Asia) but poorer agreement in others (e.g.equatorial rainforest margins, Siberian larch belt), depending on the degree of primary N limitation and the strength and direction of a shift between herbaceous and woody PFTs.
While a number of N-enabled DVMs now exist, representations of the N cycle and N-C coupling differ substantially between models in emphasis, degree of detail as well as fundamental assumptions made, reflecting different modelling goals, technical constraints, but also limitations to understanding of how the N cycle functions on a detailed, quantitative, and globally generalisable level (Zaehle and Dalmonech, 2011).The approach adopted in LPJ-GUESS is relatively detailed with respect to the N cycling and responses of vegetation.A fundamental assumption is made that plants through plastic response mechanisms acquired during evolution seek to optimise individual performance by adjusting their resource uptake and utilisation -within prescribed limits -along climatic gradients and as local conditions change (Field et al., 1992).Three such response mechanisms by plants are explicitly represented in the model: the response of stomatal conductance to CO 2 and soil water, the response of rubisco capacity to light and N availability, and the relative allocation of C and N to roots versus foliage, depending on experienced soil water and N-mediated stress.At a higher organisational level, differential availability of CO 2 , light, water and N affect the outcome of competition and rate of approach to a state of competitive exclusion between larger (taller) and smaller woody plant individuals and between co-occurring PFTs.Geographic gradients and temporal variability in C and N balance in the model are thus the emergent outcome of whole-system dynamics, modulated by the long-term evolution of vegetation structure and PFT composition, and will not necessarily reflect imposed representations of the primary physiological responses to external drivers.An example concerns the impact of N cycling on global C storage in the future-climate experiment.In general, studies based on N-enabled DVMs have concluded that any CO 2 -and climate-change-mediated increase in terrestrial ecosystem C storage will be smaller when N-C coupling is taken into account (Thornton et al., 2007;Sokolov et al., 2008;Churkina et al., 2009;Zaehle et al., 2010;Goll et al., 2012).The latter may be expected if adjustments in NPP dominate changes in whole-system C balance, and N mineralisation rates in temperature-and water-limited ecosystems strongly constrain any increase in NPP (Thornton et al., 2007;Zaehle and Dalmonech 2011).In the future-climate experiment LPJ-GUESS, by contrast, simulated a ca. 25 % greater future increase in C storage with N cycling enabled.This may be explained in part by the differential N responses of woody PFTs and grasses which lead to a gradual shift towards a larger woody vegetation component in areas where moderate drought reduces N availability, creating a transient sink for C as biomass accumulates in the stems of growing trees and shrubs.This mechanism outweighs the "primary", negative effect of N limitation on NPP which dominates the response seen in other studies.A closer analysis of the mechanisms of response of C and N balance to future-climate and CO 2 forcing is available in Wårlind et al. (2014).
A serious limitation of the current N cycle in LPJ-GUESS is the simplistic representation of biological N fixation, thought to be one of the primary processes that limit productivity of ecosystems (Luo et al., 2004), as a linear function of evapotranspiration.In addition to predicting substantially lower global values compared with observation-based estimates, a concern with the adopted parameterisation from Cleveland et al. (1999) is that it is phenomenological rather than mechanistic in nature (evapotranspiration is a largescale covariate of BNF -both tend to be higher in areas characterised by a warm climate with ample moisture -but not a direct driver) and can not therefore be assumed to hold true when the model is applied beyond the range of the observations on which it is based: for example, under a high-CO 2 future global climate (Wang and Houlton, 2009).The implementation of BNF as dependent on simulated actual evapotranspiration (AET) may also be questioned, as the relationship between AET and BNF, even within the range of modern observations of these two variables, may potentially change in the future as a result of down-regulation of stomatal conductance under elevated [CO 2 ], increasing the water-use efficiency of ecosystems, and due to differential responses of Nfixing versus non-fixing plants to elevated [CO 2 ] (Ainsworth and Long, 2005).A parsimonious alternative would be to force the model with a global climatology (time-independent map) of BNF, constructed for example from the Cleveland et al. (1999) correlation with AET, as performed in the O-CN model (Zaehle and Friend, 2010).
This study has focused on the role of soil-N-mediated constraints on plant productivity in controlling ecosystem functioning and vegetation structure.The incorporation of N cycling in LPJ-GUESS also provides the basic structure for the modelling of ecosystem production and emission of N-based trace gases (e.g.N 2 O, NO x ), allowing the model to be used as a tool in assessments of how these climatically important gases respond to climate change and changes in N input in a consistent, process-based manner (Xu-Ri and Prentice, 2008).Feedbacks related to changes in climate, CO 2 concentration, and terrestrial N 2 O emissions can be quantified in a coupled and consistent manner (Zaehle et al., 2010).Simulating soil NO x emissions within a dynamically changing environment, accounting for factors beyond soil water content or temperature (Ganzeveld et al., 2002), facilitates the conjoint analysis of terrestrial greenhouse gas emissions, and emissions of ozone and aerosol precursors in a common framework (Arneth et al., 2010).
At the aggregate global scale, accounting for N cycle dynamics did not clearly change or obviously improve the abil-ity of the model to reproduce C and N fluxes and stocks, compared with observational studies and results from other models.However, gradients in NPP along global temperature and water availability gradients were clearly changed by the incorporation of N cycle dynamics -and these may in turn be assumed to have cascading impacts on all pools and fluxes of vegetation and soils.Such similar behaviour at global scale combined with greater contrasts between regions when N cycle dynamics are accounted for reveals the presence of compensating regional errors, cancelling one another at the global scale, in the C-only or C-N version of the model (or both).In practice, it is safe to assume that all globally parameterised models, like LPJ-GUESS, have undergone some degree of tuning to attain an overall "acceptable" fit to global metrics such as NCB, GPP/NPP, vegetation and soil C (in LPJ-GUESS, such tuning has been limited to adjustment of the global quantum efficiency scalar, as discussed above and in Appendix C).As any model, by definition, is a simplification of the real-world system it represents and thus must entail some errors, such global-scale tuning will almost inevitably trade errors in different regional contexts -where different drivers and mechanisms dominate whole-system dynamics -against one another, potentially resulting in larger errors in any particular region, even while performance at the aggregate global scale appears to be improved by tuning.It is our hope that the incorporation of N cycling in LPJ-GUESS reduces regional errors by more adequately representing the key processes underlying vegetation and ecosystem dynamics in the major climate zones of the world.In general, the intercomparison of the simulation results to benchmarks of structure and function from a range of ecosystem types in the present study provides some confidence that this may be so.Table A1.Classification scheme for deriving vegetation classes (biomes) from PFT abundances for construction of Fig. 2 (modified from Hickler et al., 2006).

Vegetation class 11
Tree LAI 1 Grass LAI Appendix B

Overview of LPJ-GUESS
Here we provide an overview of the conceptual approach, assumptions and main simulated processes common to the C-N and standard C-only versions of LPJ-GUESS.Full details are available in other publications as cited below.Soil and litter C dynamics, which have been comprehensively revised in the C-N version of the model, are not included in this overview.The new CENTURY-based scheme incorporated in the C-N version of the model is presented in Appendix C.

B1 Vegetation structure and dynamics
Vegetation dynamics in LPJ-GUESS are based on the gapphase dynamics theory as conceptualised in so-called forest gap models (Bugmann 2001).The implementation in LPJ-GUESS is fully described by Smith et al. (2001), Hickler et al. (2004) and Wramneby et al. (2008).Multiple plant functional types (PFTs; Table B1) are represented and may cooccur within simulated stands or grid cells.The basic structural unit is an average individual plant, whose state is defined by foliage and fine root compartments, for woody PFTs also sapwood and heartwood (Fig. B1).Each compartment consists of a quantity of C biomass expressed on a patch ground area basis, kg C m −2 .All model description and results in the present paper refer to "cohort mode", in which, for woody PFTs, each average individual is associated with a cohort (age class) of individuals that are born in the same year and retain the same size and form as they grow.Stem density (indiv.m −2 ) is an additional state descriptor for a cohort in a patch (see below).For herbaceous PFTs ("grasses") a single average individual represents the entire population of a given patch; cohorts are not distinguished and population density is not defined.For woody PFTs, height, stem diameter and crown area of individual trees (or shrubs) can be derived from the state variables, i.e. sapwood biomass and stem density (Fig. B1), based on PFT-specific allometric constants (see Sect.B4).
Patches (0.1 ha) correspond to the maximum expected area of influence, in terms of competition for light and soil resources, of a large individual tree on its neighbours.Replicate patches are simulated to account for landscapescale heterogeneity in stand structure and composition arising from the differential disturbance histories of different patches (Fig. B1).In a standard simulation, generic patchdestroying disturbances (representing, for example, windstorms or landslides) recur stochastically with a prescribed expected frequency, typically 0.01 yr −1 .Wildfires are also simulated (Thonicke et al., 2001) and result in partial destruction of the biomass of an affected patch.
Bioclimatic limits (Sitch et al., 2003) define the environmental envelope within which PFTs are able to establish and/or survive under current climate conditions in a sim-ulated site or grid cell (Table B1).Establishment is also affected by stand structure, the density of new individuals born each year declining as crowding reduces potential forest-floor NPP below the optimum possible in the absence of standing vegetation, preempting light and soil resources.PFTs differ in terms of a suite of prescribed parameters that together define their life-history strategy.In general, shade-tolerant and shade-intolerant PFTs are distinguished (Hickler et al., 2004; Tables B1 and B2).Shadeintolerant PFTs exhibit relatively dense establishment under initial, non-crowded conditions, rapid height growth, a short average life span, and poor tolerance of resource stress (Table B2), the lattermost resulting in markedly increased mortality, reducing stem density and biomass, under shading or other resource-deficit-mediated stress.Shade-tolerant PFTs have lower initial establishment and slower height growth but a longer expected life span and lower mortality under shading or resource stress (Table B2).Interactions between PFTs differing in life history and other characteristics greatly influence the structural and compositional evolution of stands (patches) following a disturbance; under climate conditions promoting sufficient productivity for stand closure, a classic secondary successional series, with shade-intolerant trees displacing an initial layer of grass, to be subsequently succeeded by shade-tolerant trees, typically emerges (Hickler et al., 2004).

B2 Vegetation C and water balance
Exchange of CO 2 and water vapour by the vegetation canopy is computed on a daily time step by a coupled photosynthesis and stomatal conductance sub-model based on the Collatz et al. (1991Collatz et al. ( , 1992) ) simplification of the Farquhar biochemical model, with upscaling from leaf to canopy level following the strong optimality approach of Haxeltine and Prentice (1996a, b).In this model, photosynthesis, net of photorespiration, is the smaller of an electron-transport-limited and a carboxylation-limited rate (Collatz et al., 1991), and is affected by incoming photosynthetically active radiation (PAR), temperature, intercellular [CO 2 ] and carboxylation capacity, V max .The latter is determined prognostically based on the assumption (Haxeltine and Prentice, 1996a) that plants allocate N (for investment in the enzyme rubisco that determines V max ) throughout the canopy in a manner that maximises net CO 2 assimilation at the canopy level (see further discussion in Appendix C).A scaling coefficient, α a , applied multiplicatively to leaf-level daily gross photosynthesis from the biochemical model, accounts for reduction in quantum efficiency (CO 2 assimilation per unit PAR) from the leaf to canopy level.An estimate of 0.5 was suggested by Haxeltine and Prentice (1996b) and is also adopted in standard LPJ-GUESS.In addition to accounting for quantum efficiency reductions due to spectral factors such as scattering and absorption of PAR by non-photosynthetic plant surfaces, α a likely compensates for the absence of nutrient limitation in C-only Stomatal conductance (g c ) influences water vapour loss (transpiration) from the canopy and intercellular [CO 2 ], thereby coupling C and H 2 O cycling by plants and ecosystems.Aggregate stomatal conductance at the canopy scale is determined by jointly solving the biochemically based expression for photosynthesis (described above) and an alternative expression that relates photosynthesis to g c through the diffusion gradient for CO 2 implied by the ratio of intercellular to external (i.e.atmospheric) [CO 2 ] ( Haxeltine and Prentice, 1996b).
Plants are subject to maintenance and growth respiration, which are deducted from gross photosynthesis to derive NPP.Leaf respiration scales linearly with V max (Haxeltine and Prentice, 1996a).For the remaining living tissue compartments, i.e. fine roots and sapwood, maintenance respiration depends on N content and follows a modified Arrhenius-dependency on temperature (Lloyd and Taylor, 1994).Growth respiration is one-third of NPP (Ryan 1991).
Evapotranspiration (ET) encompasses transpiration by plant canopies, evaporation from exposed soil surfaces and evaporation of water intercepted by plant canopies during precipitation events.Canopy transpiration under demandlimited conditions is related to g c based on an empirical boundary layer parameterisation (Huntingford and Monteith, 1998) that expresses large-scale ET as a hyperbolic dependency on surface resistance (the inverse of g c ), thus avoiding the need for humidity as a driving variable for the model.Under supply-limited conditions, i.e. when plant water uptake from the soil is unable to meet the demand for water vapour by the atmosphere, ET is determined as a proportion of a maximum rate scaled by root zone water uptake (Haxeltine and Prentice, 1996b;Sitch et al., 2003).Soil evaporation and canopy interception are modelled as described in Gerten et al. (2004).

B3 Soil hydrology
LPJ-GUESS employs a two-layer "leaky bucket" soil hydrology scheme with percolation between layers and deep drainage.Details are available in Gerten et al. (2004).The soil of a particular grid cell may be classified into one of nine possible texture classes, affecting water holding (field) capacity, percolation and thermal properties (Sitch et al., 2003).The upper and lower soil layers are of 0.5 and 1 m depth, respectively.Incoming rainfall enters the upper soil layer, replenishing plant-available soil water up to field capacity.Excess water is exported as surface runoff.Percolation from the upper to the lower soil layer (Neilson 1995) supplies the lower layer with water.
Table B1.PFT characteristics and parameter values adopted in this study.Additional parameters specific to the C-N version of LPJ-GUESS are presented in Appendix C. Parameters common to the climate zone, shade tolerance, leaf type and growth form categories distinguished are shown in Table B2.T c,min = minimum coldest-month temperature for survival and establishment; T c,max = maximum coldest-month temperature for establishment; GDD 5 = minimum degree-day sum above 5 • C for establishment; r fire = fraction of individuals surviving fire; a leaf = leaf longevity; a ind = maximum expected non-stressed longevity.
PFT 1  Phen-Climate Shade Leaf Growth T c,min T c,max GDD 5 r fire a leaf a ind ology 2 zone 3 tolerance 4 type 5 form 6 ( 1 BNE = boreal needleleaved evergreen tree; BINE = boreal shade-intolerant needleleaved evergreen tree; BNS = boreal needleleaved summergreen tree; TeBS = temperate broadleaved summergreen tree; IBS = temperate shade-intolerant broadleaved summergreen tree; TeBE = temperate broadleaved evergreen tree; TrBE = tropical broadleaved evergreen tree; TrIBE = tropical shade-intolerant broadleaved evergreen tree; TrBR = tropical broadleaved raingreen tree; C3G = C 3 (cool) grass; C4G = C 4 (warm) grass.A simple snowpack is simulated, with precipitation adding to snow storage when the air temperature is below freezing point.Snowmelt is simulated as a linear function of positive air temperature and contributes water to the upper soil layer.
PFTs differ in their prescribed relative access to water from the upper and lower soil layers (Table B1), affecting stomatal conductance -and thereby photosynthesis and individual and population growth -in different hydrological regimes (Hickler et al., 2006).

B4 Plant growth and phenology
After deduction of 10 % for reproductive costs, plants invest their annual NPP in biomass growth.For woody PFTs, relative allocation to the biomass compartments leaves, fine roots and sapwood is determined at the level of the average individual of a cohort in a patch by satisfying mechanical, functional balance and demographic constraints.Stem diameter increases with height (Huang et al., 1992): where H is stem height (m), D is stem diameter (m) and k 2 is a PFT-specific constant (Table B2).Sapwood -which contains the transport vessels that supply leaves with water and the whole plant with assimilates of photosynthesis -is assumed to maintain a cross-sectional area in proportion to the total area of leaves (MacDowell et al., 2002): where LA is annual maximum individual leaf area (m 2 ), SA is sapwood cross-sectional area (m 2 ) and k LA:SA is a PFTspecific constant (Table B2).Relative investment in aboveand below-ground resource-uptake surfaces -leaves and fine roots -is affected by resource stress experienced the preceding growing season: where C leaf and C root denote the annual maximum C biomass of leaves and fine roots, respectively (kg C m −2 ); lr max is a PFT-specific constant ( 1 High values indicating strongly reduced establishment as growth conditions at the forest floor become unfavourable, e.g. as a result of low PAR levels (Fulton 1991;Hickler et al., 2004).
a fourth allometric constraint, linking crown area (the inverse of density assuming perfect space-filling by tree crowns) to stem diameter: where CA is crown area (m 2 ) and k 1 is a PFT-specific constant (Table B2).Allocation for herbaceous PFTs, which lack stems in the model, is subject only to the functional balance constraint (Eq B3).
A functional trade-off between leaf quality and short-term carbon return is represented as a fixed dependency of specific leaf area (SLA, leaf area per unit C mass, m 2 kg C −1 ) on leaf longevity, based on a well-documented global relationship between these variables (Reich et al., 1992(Reich et al., , 1997)).Longer-lived leaves carry a larger cost in terms of allocated C per unit leaf area, but provide a more-than-equal return in terms of assimilated C if temperatures and soil conditions are conducive to achieving positive C balance for a sufficiently long growing season each year.A single global relationship based on Reich et al. (1997) in the standard C-only version of LPJ-GUESS (Sitch et al., 2003) was replaced in the C-N version of the model by separate relationships for needleleaves and broadleaves from Reich et al. (1992).This change has consequences for C and N dynamics (see Appendix C).
Shed leaves and fine roots become part of the above-ground and soil litter pools where they enter the soil decomposition cycle (see Appendix C).Sapwood is converted to non-living heartwood, which contributes to stem dimensions but is no longer involved in water and solute transport, and does not respire.An explicit phenological cycle is simulated only for leaves in summergreen and raingreen PFTs.Summergreen phenology is based on a PFT-dependent accumulated growing degree-day (GDD) sum threshold for leaf onset, with leaf area rising from 0 to the pre-determined annual maximum (LA) linearly with an additional 200 degree days above a threshold of 5 • C. For summergreens, growing season length is fixed, all leaves being shed after the equivalent of 210 days with full leaf cover.Raingreen PFTs shed their leaves whenever the water stress scalar ω (see above) drops below a threshold, ω min , signifying the onset of a drought period or dry season.New leaves are produced, after a prescribed minimum dormancy period, when ω again rises above ω min .

B5 Plant functional type parameters
The PFTs adopted and the characteristics and parameter values distinguishing them in the simulations of the present study are shown in Tables B1 and B2 (Parton et al., 1993(Parton et al., , 2010)).
f (S) is a dimensionless scalar in the range 0-1 determined from soil fractional silt plus play content (S) following Parton et al. (1993): Litter resulting from vegetation turnover (mortality or phenology), effected in the model on the last day of a given year, is transferred to the litter SOM pools on the first day of the following year.Leaf and root litter is partitioned into structural (resistant to decomposition) and metabolic (readily decomposable) fractions based on the estimated lignin : N ratio (Parton et al., 1993): where F m and (1 − F m ) are the metabolic and structural litter fractions, respectively; λ is assumed lignin content as a fraction of total C mass (leaves: 20 %; fine roots: 16 %); and cn is the prognostic C : N ratio of the incoming material.Sapwood and heartwood biomass lost due to mortality or disturbance enter the fine and coarse woody debris litter pools, respectively.

C4 Vegetation N cycling
Plants obtain N for allocation to their biomass compartments leaves, fine roots and (for woody PFTs) sapwood through root uptake from the soil mineral N pool N avail .N uptake takes place daily and is the smaller of current supply, i.e.N avail , and demand, subject to a maximum constraint on total uptake.
Vegetation N demand is based on a solution of the carboxylation capacity of rubisco (V max ) that maximises net assimilation at the canopy level given current temperature, light interception and intercellular [CO 2 ], the latter affected by ambient [CO 2 ] but also by the influence of soil moisture and boundary layer humidity on stomatal conductance (Haxeltine and Prentice, 1996a, b;Sitch et al., 2003).Following Haxeltine and Prentice (1996a) where T is air temperature ( • C); L is day length (s); C leaf is leaf C mass, accounting for canopy phenology (g C m −2 ); and f (LAI) is a modifier dependent on current leaf area index (LAI, m 2 m −2 ) that accounts for the empirical finding that leaf N content declines more gradually with canopy depth compared to incoming sunlight (Lloyd et al., 2010;Peltoniemi et al., 2012): Based on Eq. (C10), the target leaf C : N mass ratio may be calculated as C leaf / N leaf .Leaf C : N is, however, constrained to remain within prescribed boundaries [CN leaf,min , CN leaf,max ] based on observations (Reich et al., 1992;White et al., 2000).C : N ratios for the further compartments' fine roots and (for woody PFTs) sapwood are assumed to vary proportionately with leaf C : N, fine roots maintaining a C : N ratio 1.16 times higher, and sapwood 6.9 times higher, than leaves (Friend et al., 1997).Since allocation of the current year's NPP is effected only once per year in LPJ-GUESS, allocation ratios (proportion of biomass increment allocated to each respective compartment) from the previous year are assumed when computing daily demand for allocation to fine roots and sapwood.Plants maintain a store of labile nitrogen, N store (kg N m −2 ), to buffer fluctuations in the balance between N demand and supply from the soil mineral N pool.Following Zaehle and Friend (2010), the maximum capacity of N store is related to current size as where C sap , C root , C leaf and N leaf denote sapwood C mass, fine root C mass, leaf C mass and leaf N mass, respectively, on allocation the previous year; k is set to 0.05 for evergreen woody, 0.15 for deciduous woody and 0.3 for herbaceous PFTs.The store is replenished by uptake from N avail .
Daily N demand for allocation to N store is computed as max 0, NPP a NPP y−1 (N store,max − N turnover ) − N store , (C13) where NPP a is the NPP accumulated since the beginning of the current year and NPP y−1 is the NPP of the previous year; are N turnover is the expected amount of N to be reallocated from turnover of leaves, fine roots and sapwood, based on tissue C : N ratios and biomass from the previous year.
Where the current day's bulk N demand cannot be fulfilled by the present size of N avail , total uptake is reduced to N avail .Uptake may also be further reduced to a ceiling, N up,max , computed following Zaehle and Friend (2010) as N up,max = 2N up,root f (N avail )f (T soil )f (NC plant )C root , (C14) where N up,root is a linear scalar of the maximum N uptake per unit fine root biomass, C root , assuming a proportional increase in uptake capacity with root exploration volume, assigned the fixed values 2.8 and 5.51 g N kg C −1 day −1 for woody PFTs and grasses, respectively (Rothstein et al., 2000;Macduff et al., 2002).Modifiers account for the effects of the current mineral N pool, soil temperature (Eq.C5) and plant N status on uptake capacity, as follows: representing a combined linear and saturating effect of mineral N concentration on N uptake (Zaehle and Friend, 2010), with k m , the half-saturation concentration for N uptake, set to 1.48 g N m −3 for woody PFTs (Rothstein et al., 2000) and 1.19 g N m −3 for grasses (Macduff et al., 2002); z soil is the soil column depth (1.5 m).Vegetation N demand and uptake are computed each daily time step for each average individual plant (in practice, each age/size class in each replicate patch for woody PFTs, and once for the herbaceous ground layer in each patch).In the event that bulk demand cannot be met by the available N supply, the supply is partitioned among individuals in proportion to their relative uptake strength, f Nup , which is related to estimated fine root surface area following where X (indiv.m −2 ) is stem (cohort) density (included as a weighting factor for the most abundant cohorts) and k Nup , set to 1.6 for woody PFTs and 1.9 for grasses, weighs N uptake towards PFTs having shallower root distributions, coinciding with an assumed greater concentration of available N in the upper soil layer (Franzluebbers and Stuedemann, 2009).Equation (C18) also implies that plants become more efficient at taking up N when their store of relatively mobile N approaches its lower limit, e.g. through physiological up-regulation of root uptake capacity (Raynaud and Leadley, 2004); the existence of such a response is also suggested by studies of plant-mycorrhizal associations, which are often more developed in N-depleted habitats (Olsrud et al., 2004).
Where N uptake is insufficient to meet individual demand, individuals attempt to fulfil the deficit using their current labile N store.If demand is still not met after the N store is depleted, rubisco capacity and thereby leaf and whole-plant demand (as well as photosynthesis) are reduced to the maximum level that can be satisfied given the current supply plus storage, effecting N limitation.The N store is replenished, up to its maximum capacity (see above), on the last day of the year by retranslocation of up to 50 % of the N mass of shed leaves, fine roots and sapwood on conversion to heartwood (Aerts 1996, Vergutz et al., 2012).

C5 Plant growth and C and N allocation
Plant growth takes place on the last day of the simulation year by allocation of annual accrued NPP to the biomass compartments leaves, fine roots and (for woody PFTs) sapwood subject to allometric constraints.The allocation scheme in the standard C-only version of LPJ-GUESS is described in Appendix B. The only modification resulting from the incorporation of N cycling in the model is the addition of an N stress scalar (υ) in the functional balance constraint that governs the relative allocation of biomass to foliage versus fine roots: where lr max is a PFT-specific constant (Table B2) and ω is a soil moisture stress scalar in the range 0-1, with smaller values reflecting increased soil moisture stress (Sitch et al., 2003;Appendix B).Where N stress exceeds soil moisture stress, this results in an increased allocation of biomass to fine roots at the expense of foliage: where CN leaf,aopt is the leaf C : N ratio that would have been realised if plant N demand had been fulfilled by the available supply plus storage every day of the current year; C leaf and N leaf are realised C and N mass, accounting for N limitation.

C6 Additional updates
The incorporation of N limitation naturally results in a reduction in simulated NPP, regionally and globally, relative to the C-only version of the model, which lacks such limitation.To compensate for this nutrient effect on global C balance and fluxes, the quantum efficiency scalar α a (see Appendix B) was recalibrated to a value that resulted in simulated global C fluxes within the approximate range of observation-based estimates (Table 1).For the C-only simulations of this study, which were performed with the C-N model, but with N limitation switched off, α a was likewise calibrated to the global fluxes.The resulting settings of this parameter were 0.70 and 0.55 with N limitation enabled and disabled, respectively.It may be postulated that the 15-percentage-point differential between these values corresponds to the global limitation of primary production attributable to N limitation, whereas the residual difference of 30 % between realised and potential canopy quantum efficiency, obtained with N limitation enabled, more closely reflects the spectral factors traditionally invoked to explain this parameter.
In conjunction with the incorporation of N cycling, an equation linking SLA to the PFT parameter leaf longevity (a leaf ; Table B1), originally adopted from Reich et al. (1997), was replaced with separate parameterisations of the same relationship for needleleaved and broadleaved PFTs, following Reich et al. (1992).The new equation has the form SLA = 0.2g(β 0 + β 1 log 10 12a leaf ) g(p) = 10 p , (C21) with SLA in m 2 kg C −1 and a leaf in yr.The regression coefficients {β 0 , β 1 }, fitted to a global data set by Reich et al. (1992), are set to {2.41, −0.38} and {2.29, −0.40} for needleleaves and broadleaves, respectively.As a consequence of this update, leaves are simulated to be generally thicker, with lower SLA and consequently reduced PAR per unit invested leaf C. The global data presented in both the 1992 and 1997 Reich et al. (1992) papers are, however, more faithfully reproduced, suggesting the presence of a unit conversion error in the original implementation.The resulting reduction in productivity per unit leaf C more strongly penalises species with short-lived leaves, particularly deciduous species and grasses, providing one explanation for an increased dominance by woody PFTs relative to grasses in simulations with the updated model, whether or not N limitations are enabled (Sect.3.2).

Fig. 1 .
Fig. 1. (A) Mean annual NPP (kg C m −2 yr −1 ) for 1961-1990 simulated by the C-N version of LPJ-GUESS in the global hindcast experiment; (B) N limitation of NPP, showing the value obtained by dividing mean NPP in each grid cell for 1961-1990 from the C-N simulation by the value obtained in a C-only simulation with the same value of the global quantum efficiency scalar α a (see text).Values close to 1 thus denote "no limitation", lower values increasing limitation of NPP by N availability.

Fig. 2 .
Fig. 2. Potential natural (Haxeltine and Prentice, 1996b; Hickler et al., 2006) (A) and simulated vegetation for 1961-1990 from the (B) C-N and (C) C-only simulations.The classification scheme for transforming PFT abundances from the model simulation to the vegetation classes (biomes) in the legend modified from Hickler et al. (2006) (see Appendix A, TableA1).

Fig. 3 .
Fig. 3. Comparison of data on forest structure (A: tree density; B: height), (C) NPP and (D) N use efficiency (NPP per unit plant N uptake) for the CANIF European forest sites with simulation results from LPJ-GUESS.Model runs were adjusted to match the observed tree density at each site in (A).All other quantities are simulated without site-specific calibration.(C) includes biome-average data for 132 sites from the Luyssaert et al. (2007) database.

Fig. 4 .
Fig. 4. Size structure and generalised composition of boreal forest stands along a chronosequence constructed by ranking 12 virgin forest stands in northern Sweden in order of age.(a) Observed data modified from Linder et al. (1997); (b) average for corresponding locations simulated under repeated 1961-1990 climate.Shade-intolerant = aggregated data for Betula spp., other deciduous, and Pinus sylvestris; shade-tolerant = Picea abies.Note: x axis not to scale for (a); each tick mark corresponds to a specific stand arranged in rank order of age along the axis.
hancement from high to middle latitudes is steeper in the new C-N version of the model.This is explained by N limitation at the cold end of the gradient, where the N-mineralisationand -fixation-limited supply of soil N is unable to match an increased demand from the less [CO 2 ]-limited canopy.NPP enhancement drops to < 10 % over much of the Boreal Zone and Arctic, compared with 5-15 % in the C-only model.Similar findings were reported from a comparison between the C-N and C-only versions of the O-CN model by Zaehle et al. (2010).

Fig. 5 .
Fig. 5. Component biomass allometry for broadleaf-dominated (A, C) and needleleaf-dominated (B, D) grid cells from C-N and C-only simulations superimposed on data for forest stands from the Cannell (1982) and Usoltsev (2001) databases.Curves are fitted to model data from the 0.5 • × 0.5 • grid cells encompassing stand coordinates from the observed databases.Allometric relationships from the land surface models considered by Wolf et al. (2011) (grey curves) are reproduced from that paper.N = tree population density (ha −1 ); M = total biomass of average tree (kg); M fol = foliage biomass of average tree (kg); M stem = stem biomass of average tree (kg).

Fig. 6 .
Fig. 6.Geographic pattern of simulated NPP enhancement (%) (1996-2002) resulting from a step increase of atmospheric CO 2 concentration from ambient to 550 ppmv; (A) pattern of NPP enhancement from C-only model in Hickler et al. (2008); (B) corresponding pattern from LPJ-GUESS C-N (present study); (C) average NPP enhancement by latitude band from Hickler et al. (2008) and present study.

Fig. 7 .
Fig. 7. Nitrogen required to support terrestrial carbon uptake between 2000 and 2100 simulated in the C-N and C-only simulations under CO 2 -only and CO 2 -and-climate-change scenarios based on RCP 8.5 radiative forcing of the MPI-ESM-LR general circulation model (GCM).High and low N supply limits proposed by Hungate et al. (2003) are shown.For C-only the upper N-requirement value assumes a fixed tree C : N of 200 following Hungate et al. (2003), the lower value that all new tree carbon is allocated to wood with a C : N of 500.Soil is assumed to have a fixed C : N of 15.For C-N, N requirements are as simulated by the model.
f (NC plant )= max 0, NC plant −1/CN leaf,min 2/(CN leaf,max+CN leaf,min )−1/CN leaf,min ,(C16) representing a tendency for N uptake to increase as the concentration of relatively mobile N compounds within the plant, characterised by NC plant (below), declines.CN leaf,min and CN leaf,max are the prescribed minimum and maximum bounds for leaf C : N. NC plant = N leaf + N root C leaf + C root (C17) et al.: Implications of incorporating N cycling and N limitations in an DVM

Table 1 .
Mean global C and N stocks and fluxes simulated for 1961-1990in the global hindcast experiment, and according to literature estimates.

Table 2 .
Slope (β) of the scaling relationships depicted in Fig. 5 for observed forest stands and for the C-N and C-only simulations.For symbols see caption of Fig. 5.

2027-2054, 2014 2050 B. Smith et al.: Implications of incorporating N cycling and N limitations in an DVMTable B2 .
PFT characteristics common to the shade tolerance, leaf type, growth form and climate zone categories in TableB1.