Accounting for Carbon and Nitrogen interactions in the Global Terrestrial Ecosystem Model ORCHIDEE (trunk version, rev 4999): multi-scale evaluation of gross primary production

. Nitrogen is an essential element controlling ecosystem carbon (C) productivity and its response to climate change and atmospheric [CO 2 ] increase. This study presents the evaluation – focussing on gross primary production (GPP) – of a 15 new version of the ORCHIDEE model that gathers the representation of the nitrogen cycle and of its interactions with the carbon cycle from the OCN model and the most recent developments from the ORCHIDEE trunk version. We quantify the model skills at 78 Fluxnet sites by simulating the observed mean seasonal cycle, daily mean flux variations, and annual mean average GPP flux for grasslands and forests. Accounting for carbon-nitrogen interactions does not substantially change the main skills of ORCHIDEE, except for the site-to-site annual mean GPP variations, for which the 20 version with carbon-nitrogen interactions is in better agreement to observations. However, the simulated GPP response to idealized [CO 2 ] enrichment simulations is highly sensitive to whether or not carbon-nitrogen interactions are accounted for. Doubling the atmospheric 2 ] rise, land use change (LUC) and evolution of nitrogen atmospheric deposition (Ndep), nitrogen synthetic fertiliser (Nfert) and manure (Nman) on a yearly basis. CRU-NCEP data for the period was used for simulating the period from 1861 to 1901, while from 1901 onwards the climatic forcing from the matching year was used. This simulation, which accounts for carbon-nitrogen interactions, is named S1-CNdyn.

carbon-nitrogen interactions are not accounted for, and that carbon-nitrogen interactions are essential in understanding global terrestrial ecosystem productivity.

Introduction
Global Terrestrial Ecosystem Models (GTEMs) are mathematical models which are dedicated to provide a better understanding of terrestrial ecosystem functioning and its interplay with environmental drivers such as temperature or 5 precipitation. GTEMs aim at simulating the spatial patterns of the fluxes of carbon, water and energy between the land surface and the atmosphere, and their time evolution, in particular in a context of climate change. For over a decade, GTEMs account for climate forcing as well as the effect of atmospheric CO 2 concentration (atmospheric [CO 2 ]) on ecosystem productivity (Melillo et al., 1995). Atmospheric [CO 2 ] is a key driver of the assimilation of carbon by photosynthesis. While the current atmospheric [CO 2 ] value is nearby the optimal value for carbon assimilated by plants with a C4 photosynthetic 10 pathway, it is still suboptimal for plants with a C3 pathway . In the context of global change where atmospheric [CO 2 ] is increasing, quantifying the so-called [CO 2 ] fertilization effect, i.e., the increase in ecosystem productivity associated to increasing atmospheric [CO 2 ] has been on the forefront . 15 Most of the GTEMs estimate a large global land carbon sink over the 21 st century when used in Earth System Models for climate predictions (of the order of 120-270 GtC over the period 2010-2100, depending on the Representative Concentration Pathways), mainly driven by atmospheric [CO 2 ] increase . Even if atmospheric [CO 2 ] will be plentiful in the future, it remains questionable whether sufficient nutrients, in particular nitrogen, will be available to fully sustain the increase of primary production associated solely to the rise of [CO 2 ]. A recent study  estimated that 40 to 20 80% of the carbon sequestration on land projected by simulations without nutrient limitations for the period 1851-2100, would not occur if nitrogen limitation and carbon-nitrogen interactions were accounted for in GTEMs embedded into Earth system models. The 40% variation in the projected N-limitation of the land carbon sink was reported to depend on the evolution of the anthropogenic production of reactive nitrogen and associated atmospheric nitrogen deposition, which differ for each Representative Concentration Pathway . 25 Only two GTEMs involved in the last exercise of the Climate Modelling Intercomparison Project (CMIP5) accounted for carbon-nitrogen interactions (CESM1-BGC and NorESM1-M). Since then many GTEMs have been further developed to account for the impact of the nitrogen cycle Esser et al., 2011;Fisher et al., 2010;Goll et al., 2017;Smith et al., 2014;Sokolov et al., 2008;Wang et al., 2007Wang et al., , 2010Xu-Ri and Prentice, 30 2008;, some of which being included in an Earth system model. Among the GTEMs that developed carbon-nitrogen interactions was an ORCHIDEE-derived model named OCN , that has been used in several studies . However, this pioneering development (2007)(2008)(2009)(2010) has not been embedded in subsequent versions of the ORCHIDEE model, especially with respect to the coupling to the French IPSL (Institut Pierre Simon Laplace) Earth system model. The present paper presents and evaluates a recent modelling effort consisting of the merge of the ORCHIDEE trunk version (r3977) with the carbon-nitrogen interactions based on . It describes all changes made in the original nitrogen cycle and carbon-nitrogen interactions, linked to major 5 updates in the water, carbon and energy budgets in ORCHIDEE since the first development of OCN, together with a thorough evaluation of simulated gross carbon uptake by plants.

Methods
Following a description of the model developments, the evaluation is focused on the carbon cycle and on the added value of including the nitrogen cycle for the purpose of simulating gross carbon uptake fluxes, including also the impact on the 10 related plant transpiration. The evaluation consists of: (1) site-level simulations in order to assess the overall performance of the ORCHIDEE model at simulating GPP flux at Fluxnet stations (annual mean value, seasonal variations, site-to-site differences); (2) sensitivity tests to quantify the contributions of accounting for seasonal and site-to-site variations of the leaf C/N ratio to the simulated seasonal variations and mean annual GPP; (3) idealized simulations to quantify the impact of Nlimitation on GPP under [CO 2 ] enrichment scenarios; and (4) global simulations in order to evaluate and analyse seasonal, 15 long-term variations and global distribution of the simulated GPP. While the OCN model , the predecessor of ORCHIDEE r4999, has already been evaluated over a restricted set of sites for which C and N data are available, the extended C flux dataset from the Fluxnet network has so far not been used for an in-depth evaluation of an ORCHIDEE version that includes the N cycle and the C/N interactions.

Model description 20
This section focuses on the major modifications that were included in the ORCHIDEE model since the first implementation of a nitrogen cycle in the OCN model ). The ORCHIDEE model calculates the exchange of energy, water, carbon and nitrogen at the atmosphere-surface interface and within the soil-plant continuum. The main modelling structure originates from Ducoudré et al., (1993) and Krinner et al., (2005) for water-and energy-related and carbon-related processes, respectively. The spatial discretization depends on the modelled process. The energy budget is computed at the 25 grid cell level, without accounting for differences between tiles within the grid cell. Water budgets are now calculated for three tiles per grid cell, one for the bare soil, one for tree cover and one for herbaceous cover. The carbon budget and related fluxes are computed for each vegetated tile within a grid cell. In ORCHIDEE, the carbon in the vegetation and soil is split in thirteen classes, based on the concept of Plant Functional Type (PFT). Different species, which share similar characteristics regarding their architecture, phenology or location, are gathered in a single PFT. The thirteen PFTs used in ORCHIDEE are 30 listed in Table 1. Compared to the model version that was used as the basis for OCN, a large portion of the code has been revised and new developments added (Peylin et al., in prep.). The main changes consist of : i. a multi-layer soil hydrological scheme which accounts for water diffusion and deep drainage, based on the initial work of de  that was shown to better simulate soil water dynamics compared to the initial double buckets scheme as described in Ducoudré et al. (1993). This feature within the current version of ORCHIDEE has 5 been recently evaluated over Amazonia  and Africa (Traore et al., 2014); ii. a revised thermodynamic scheme which accounts for the heat transported by liquid water into the soil, in addition to improvements in the representation of heat conduction process (Wang et al., 2016) which resolves former inconsistencies between the soil water and energy dynamics; iii. an analytical solution to the set of equations governing the soil organic matter (SOM) pools and their time evolution 10 driven by the litter input and the climate conditions in terms of soil temperature and humidity (Lardy et al., 2011).
The latter is now used to determine SOM pools associated to initial conditions and guarantees steady state conditions better than iterative simplified schemes (Krinner et al., 2005).

iv.
A revised parameterisation of the vegetation and snow albedo (i.e. optimized parameters using remote sensing albedo data from MODIS sensor); 15 The developments in ORCHIDEE r4999 that relate to the nitrogen dynamics within the soil-plant-atmosphere continuum and the dependency between carbon and nitrogen cycles, as well as the allocation scheme and the short-and long-term storage pools dynamic mostly follow the work of  and . The nitrogen cycle is added at the PFT level as for the carbon cycle and for each carbon pool there is a corresponding nitrogen pool, with C/N ratios 20 evolving through time. Leaf C/N ratio is dynamic and varies as a result of the nitrogen supply by roots and demand for biomass allocation (see Sect. 2.1.3 for details). C/N stoichiometry of the other living biomass pools (i.e., below and aboveground sapwood, below and above-ground heartwood, fruit and fine roots) is driven by the C/N ratio of the leaves, but multiplied by a pool-dependent factor f cn (f cn equals 1.2 for fine roots and fruit, and 11.5 for the other pools). SOM decomposition follows the scheme of  in which C/N ratios of SOM pools are expressed as a function of 25 soil mineral nitrogen content (ammonium and nitrate). This scheme is rather simple and does not account for known processes such as the priming effect (Fontaine et al., 2007). As a consequence, carbon decomposition rates are independent of the C/N ratio of the SOM pools, which facilitate the use of an analytical solution for quantifying the carbon content of SOM pools at equilibrium.

30
In ORCHIDEE r4999, the fate of mineral nitrogen in the soil follows the formalism of the OCN model ) mainly based on the DNDC model . The formalism accounts for ammonium (NH 3 /NH 4 + ), nitrate (NO 3 -), nitrogen oxides (NO x ) and nitrous oxide (N 2 O) soil pools and associated emissions due to nitrification (the oxidation of NH 3 /NH 4 + in NO 3 -) and denitrification (the reduction of NO 3 up to the production of N 2 ) processes. NO 3 -(resp. NH 4 + ) uptake by roots is modelled as a function of the NO 3 -(resp. NH 4 + ) available in the soil  and of the root biomass. The more root biomass or the higher the soil NO 3 -(resp. NH 4 + ) pool, the higher the NO 3 -(resp. NH 4 + ) uptake is. Nitrogen inputs in the soil/plant system are related to (i) atmospheric nitrogen deposition under the form of NH x and NO y components, (ii) biological nitrogen fixation (BNF) on any land category and (iii) nitrogen fertilisation over managed grasslands and croplands. The nitrogen output fluxes are associated to 5 runoff and leaching, and emissions of NH 3 , NO x , N 2 O and N 2 .
The main modifications compared to the work of  are related to the carbon assimilation or photosynthesis scheme ) and refinements of the N-dependency of the photosynthetic activity ). These developments are described in the following sections 2.1.1 and 2.1.2. Furthermore, in the present study, the 10 BNF rates are computed as a function of evapotranspiration following the approach of Cleveland et al. (1999). For this purpose, a single climatology of evapotranspiration, based on a global ORCHIDEE simulation for present-day conditions, is used in all simulations performed in this study. As a consequence, the differences in modelled GPP by the different model configurations (see below) cannot be attributed to changes on BNF, an approach we consider reasonable due to the large uncertainties associated with the estimates of BNF (Zheng et al., 2019). The forcings used for the other nitrogen input fields 15 are detailed in Sect. 2.3.4 and 2.3.5.

Carbon assimilation scheme
The updated carbon assimilation scheme used in ORCHIDEE r4999 has been proposed by . This scheme is based on the model developed by Farquhar, von Caemmerer and Berry in the FvCB model, (Farquhar et al., 1980) which predicts carbon assimilation for C3 plants as the minimum of the Rubisco-limited rate of CO 2 assimilation (Ac) and 20 the electron transport-limited rate of CO 2 assimilation (Aj). ) propose a C4-equivalent version of the FvCB model and analytical solutions to the set of equations which link the net assimilation rate, the stomatal conductance, and the intercellular CO 2 partial pressure.
ORCHIDEE r4999 retained most of the parameter values of the FvCB model as proposed by , except 25 the parameters which define the maximum rate of Rubisco activity-limited carboxylation (Vc max , µmol CO 2 m -2 [leaf] s -1 ) and the maximum rate of electron transport (e -) transport under saturated light (J max , µmol e-m -2 [leaf] s -1 ) for C3 plants, which were replaced by the formulation and parameterization proposed by . In ORCHIDEE (r4999), Vc max and J max at temperature T, are defined as: where k is either Vc max or J max , k ref is the parameter value at the reference temperature (T ref is set to 25°C, expressed in Kelvin in Eq. (1)), T l is the leaf temperature (°K), R is the universal gas constant (8.314 J K -1 mol -1 ), ΔS k an entropy factor (J K -1 mol -1 ), and H a,k and H d,k , respectively an energy of activation and deactivation (J mol -1 ). Based on the reanalysis of the temperature dependency of Vc max and J max performed by , ΔSk with k being either Vc max or J max acclimate to temperature and consequently are expressed as linear functions of the monthly mean leaf temperature (t growth , 5 °C). The ratio of Vc max,ref to J max,ref also acclimates to temperature, and  proposed to define it as a linear function of t growth . Consequently, J max,ref can be expressed as : where a rJ,V and b rJ,V are fitted parameters of the relationship between observation-based values of J max,ref /Vc max,ref and t growth , and equal to 2.59 [-] and -0.035 [°C -1 ], respectively. 10

Nitrogen dependency of photosynthesis activity
ORCHIDEE r4999 accounts for the nitrogen limitation on photosynthetic activity in a different manner than in OCN (Friend and Kiang, 2005;, by making Vc max,ref a function of the leaf nitrogen content (N l , g [N] m -2 [leaf] ) as proposed and parameterized by : where NUE ref is the nitrogen use efficiency (µmol CO 2 g -1 [N] s -1 ). NUE ref has been widely measured (Ellsworth et al., 2004;Medlyn and Jarvis, 1999;Woodward et al., 1995) and was reported to be PFT-dependent ) (see Table 1).
Following observations of vertical N-profiles in tree canopies (Thornton and Zimmermann, 2007), N l is exponentially decreasing from the top to bottom of canopy and its value at a cumulative leaf area index (L;m 2 [leaf] m -2 [ground] , starting from 20 the top of the canopy) is defined following Dewar et al. (2012), with a specific extinction coefficient (k N , value of 0.15) that differs from the one used to calculate the light profile within the canopy (value of 0.5): where N l (0) the value of N l at top of canopy is expressed as a function of N tot (g [N] m -2 [ground] ), the total canopy nitrogen 25 content, and L tot the LAI of the total canopy : As in Dewar et al. (2012), it is assumed that the variation of the leaf nitrogen content (N l ) through the canopy is due to variation of the specific leaf area (SLA, defined as the leaf area divided by the leaf mass , m 2 [leaf] g -1 [C] ), the leaf nitrogen 30 concentration (N conc ; g [N] g -1 [C] ) being constant through the canopy. It is also assumed that the SLA at the bottom of the canopy (sla fix ) is fixed. This implies that the mean SLA of the canopy (SLA mean ) is no longer a fixed value as was formerly the case in ORCHIDEE, but varies with the total leaf mass (m leaf , g [C] m -2 [ground] ). SLA mean can be written as: and L tot as: 5

Nitrogen-related model configurations
Two model configurations were developed to allow a straightforward analysis of the effect of the nitrogen cycle on plant productivity: one with prescribed leaf nitrogen concentrations and the other with leaf nitrogen concentrations varying according to nitrogen availability. In the first configuration named "CNfix", the leaf C/N ratio is fixed to a value within a prescribed range ([CN leaf,min ; CN leaf,max ], see Table 1). In this configuration the limitation of ecosystem productivity by 10 nitrogen availability is reflected by the imposed leaf C/N ratio, which is fixed for the entire simulation (see Sect. 2.4.3 for the different simulations performed with the "CNfix" configuration). In the "CNfix" configuration, the mass balance of the N cycle within the soil-plant continuum is closed by taking the nitrogen that is needed for maintaining the imposed leaf C/N ratio from the atmosphere. Rather than implying the absence of nitrogen limitation, the "CNfix" configuration implies a fixed nitrogen limitation, which will not change over time depending on the nitrogen availability. In the other configuration 15 (labelled "CNdyn"), the leaf C/N ratio is not fixed but dynamic. The variation of the leaf C/N ratio (CN leaf , g [C] g [N] -1 ) is the outcome of the N-supply from the roots vs. the N-demand to convert the assimilated carbon into leaf, wood, root and fruit tissue each with its own C/N ratio.
Irrespective of the configuration the model first calculates the nitrogen required (G Ninit , g [N] m -2 [ground] d -1 ) for satisfying the 20 new carbon allocated to the different reservoirs, G C (g [C] m -2 [ground] d -1 ) under the assumption that CN leaf does not vary .
where f i are the fractions (unitless) of carbon allocated to leaf (l), roots (r), fruit (f) and sapwood or stalks (s) and CN i are the 25 C/N ratios (unitless) for the different biomass pools at the previous time step. Assuming that the differences in C/N ratio among the different pools are fixed, they can all be expressed as functions of the C/N ratio of the leaves and the nitrogen required can be further expressed as:

30
where λ i are unitless coefficients which reflect the differences in C/N ratio of roots, fruit, and sapwood or stalks, respectively to CN leaf . Given that all the available nitrogen is stored in the labile pool (N lab , g [N] m -2 ), the model will then check whether there is sufficient nitrogen in the labile pool to satisfy the demand (if G Ninit < N lab ) resulting in a decrease of the leaf C/N ratio of the newly allocated biomass (CN leaf,alloc ) or contrary, that the labile pool does not contain sufficient nitrogen to satisfy the demand (if G Ninit > N lab ) resulting in an increase of CN leaf,alloc . For this purpose, a dynamic elasticity variable (D leaf ) is used to dampen the leaf C/N variations and to ensure that CN leaf,alloc remains within the prescribed range of variation, [CN leaf,min ; CN leaf,max ]. At each time step, CN leaf,alloc is defined from the value of CN leaf , as: 5 where The elasticity variable D max , which avoids rapid changes in CN leaf,alloc , was slightly modified compared to its initial 10 implementation in OCN  in order to have a value of 1 for D max when CN leaf equals CN leaf,min . D max is now calculated as: where NC leaf , NC leaf,max and NC leaf,min correspond to 1/CN leaf , 1/CN leaf,min and 1/CN leaf,max , respectively and ! !"# and ! !"# are 15 two empirical parameters set to 1.6 and 4.1, respectively, in order to best fit the original function (eq. 21 of the Supplementary Material of ). In the extreme case where G Ninit is greater than N lab while CN leaf equals CN leaf,max , the new carbon allocated G C is lowered in order to maintain CN leaf at the level of CN leaf,max .
When N lab is lower than G Ninit , the closer CN leaf is to CN leaf,min , the higher the nitrogen content reduction of the newly allocated biomass. On the opposite, when N lab is higher than GN init , the closer CN leaf is to CN leaf,max , the higher the nitrogen 20 content enrichment of the newly allocated biomass.

Fluxnet GPP product
The free fair-use Fluxnet LaThuile collection (http://fluxnet.fluxdata.org/data/la-thuile-dataset/) was used to evaluate the model performance at individual sites. We selected the observation-derived GPP flux based on the NEE partitioning method 25 of Reichstein et al., (2005). From the 153 sites contained in the collection, we selected sites with vegetation belonging to a single PFT, and thus excluded sites covered by a mixture of vegetation (such as savannahs or opened forests). This is to avoid having to set fractions of PFT, that are usually uncertain, which would introduce substantial uncertainty in the evaluation process. Furthermore, sites were excluded if their vegetation cover was not explicitly represented in ORCHIDEE, such as shrublands or wetlands. In addition, because the nitrogen fertiliser inputs are not reported in the Fluxnet database for managed agro-ecosystems such as grasslands and croplands, which are fertilized with manure or synthetic forms, they were also excluded from the analyses. These two filtering criteria reduced the validation set to 86 sites. Last, we removed eight sites for which the annual mean precipitation derived from in-situ measurements was highly different of the climatological mean provided by the site principal investigator and of the value derived from the ERA-interim reanalysis (see Vuichard and 5 Papale, 2015). The selected 78 sites (Table A1)  Needleleaved Forest (BoENF) and 11 C3 natural grasslands (GRAc3).

Global MTE-GPP product 10
The observation-based MTE-GPP product (Jung et al., 2011)  and that neither atmospheric [CO 2 ] nor nitrogen availability are predictors in the training, MTE-GPP has no temporal trend which may be attributed indirectly or directly to these two driving variables.

Meteorological data
In-situ meteorology, which is typically observed at individual Fluxnet stations, was used to drive the ORCHIDEE simulations at the Fluxnet stations. In-situ data were gap-filled as described in Vuichard and Papale (2015) in order to provide continuous half-hourly records of temperature, precipitation, short-and long-wave incoming radiation, wind speed and specific humidity to the ORCHIDEE model. The mean length of the meteorological data was 5 years and ranged 25 between 1 to 16 years.
Global scale simulations were driven by CRU-NCEP meteorological data, available for the 1901-2016 period. CRU-NCEP consists of 6-hourly meteorological fields from the NCEP/NCAR reanalysis at 0.5 x 0.5 degree resolution that was bias corrected with monthly CRU data. 30

Vegetation-related information
For the site simulations at Fluxnet stations, the selection of the Plant Functional Type used to characterize vegetation at each site was based on in-situ information gathered within the Fluxnet dataset, using the International Geosphere-Biosphere Programme (IGBP) classification. For the global scale simulations, the PFT distribution within each grid cell over the time period 1860-2016 was derived from History Database of the Global Environment (HYDE v3.2; Goldewijk et al., 2017) for 5 the crop and pasture extend and from  for the specification of the forest types.

Soil data
The soil-related data used for driving ORCHIDEE are texture class, pH and bulk density at 0.5 x 0.5 degree resolution.
Texture class was coming from

Nitrogen deposition data
Monthly atmospheric nitrogen deposition (NH x and NO y ) during 1860 -2014 were taken from the IGAC/SPARC Chemistry-Climate Model Initiative (CCMI, Eyring et al., 2013). Nitrogen deposition fields are available globally at a resolution of 0. 5 15 x 0.5 degree. In the CCMI models, nitrogen emissions from natural biogenic sources, lightning, anthropogenic sources, biomass burning are accounted for, as well as the atmospheric transport of nitrogen gases. The CCMI product is used in the global N 2 O Model Intercomparison Project (NMIP, Tian et al., 2018) and is the official product for CMIP6 models without interactive chemistry components. ORCHIDEE reads total (wet and dry) nitrate (NO y ) and total ammonium (NH x ) atmospheric deposition rates from the CCMI product and uses this information to drive the nitrogen cycle within the soil-20 plant continuum. In ORCHIDEE, nitrate and ammonium are added to the respective soil mineral nitrogen pools at each model's time step, omitting nitrogen interception by the canopy.

Nitrogen fertilizer data
Global simulations were also driven by nitrogen application under the form of synthetic fertilisation or manure. Synthetic nitrogen fertilizer gridded annual data from 1960 to present-day were developed within the N 2 O Model Intercomparison 25 Project  and is based on national-level data from the International Fertilizer industry Association (IFA) and the FAO. Nitrogen fertilizer application rate between 1860 and 1960 were extrapolated assuming that gridded values for 1960 were linearly reduced to the zero in 1900. For nitrogen fertilization through manure application, gridded annual data were compiled and downscaled based on country-level livestock population data from FAO . Note that the carbon and nitrogen contained in manure represents a lateral flux from grasslands to croplands. When closing the carbon 30 and nitrogen budgets or calculating the net biome production -not applicable to this study, manure and synthetic nitrogen fertilizer should thus be accounted for in different ways.

Spin-up procedure
Each simulation requires a spin-up during which the model state variables (e.g., soil and biomass carbon and nitrogen pools) 5 are put at equilibrium. Given that the time period needed to reach equilibrium by far exceeds the length of the available insitu meteorological records, the spin-up at Fluxnet stations cycles several times over the entire record of in-situ meteorological observations. The spin-up started with a 500-year long run that makes use of the semi-analytical spin-up (see Sect. 2.1) in order to get the fluxes of litter fall but also the soil organic carbon pools at equilibrium, for [CO 2 ] atmospheric concentration and nitrogen deposition of the year 1860. Following this steady-state simulation, the spin-up continued with a 10 transient simulation from 1860 up to the start of the observation period for each site, varying [CO 2 ] atmospheric concentration and the nitrogen deposition with historical data and still cycling the meteorological in-situ data.
Likewise, a semi-analytical spin-up for the global simulations was performed for the conditions of the year 1860, by using the nitrogen input data (deposition and fertiliser fields), the [CO 2 ] value and land-use of the year 1860. Because CRU-NCEP 15 is only available from 1901, data for the period 1901-1920 was used for the spin-up. As for site simulations, we performed a transient simulation varying the different fields driving the GPP evolution (see below).

Reference simulations
From the end of the transient phase at Fluxnet stations, a set of simulations (one at each Fluxnet station) was performed over the observation period with [CO 2 ] and nitrogen deposition level (from the CCMI monthly dataset) corresponding to this 20 period. These simulations for present-day conditions (pd), in which carbon and nitrogen cycles are fully coupled (CNdyn configuration), are named pd-CNdyn.
Likewise, from the global steady state corresponding to the end of the spin-up simulation, a transient simulation from 1861 to 2016 is ran, accounting for climate change, [CO 2 ] rise, land use change (LUC) and evolution of nitrogen atmospheric 25 deposition (Ndep), nitrogen synthetic fertiliser (Nfert) and manure (Nman) on a yearly basis. CRU-NCEP data for the period 1901-1920 was used for simulating the period from 1861 to 1901, while from 1901 onwards the climatic forcing from the matching year was used. This simulation, which accounts for carbon-nitrogen interactions, is named S1-CNdyn.

Sensitivity test simulations
At each Fluxnet station, different simulations were performed in order to test the sensitivity of ORCHIDEE to different processes or forcing variables. As described in Eq. (3), leaf nitrogen content and leaf C/N ratio affect the carbon assimilated by photosynthesis, which in turn affects the leaf area index and, feeds back to carbon assimilation. Consequently, a set of two simulations was ran to investigate the impact of constraining the leaf C/N ratio on GPP, respectively in time and in time 5 and space. One simulation is named pd-CNfix-time (based on the CNfix configuration) in which the leaf C/N ratio was fixed to the time-average value at site level, which was in turn inferred from the pd-CNdyn reference simulation. In the other simulation named pd-CNfix-timePFT simulation (also based on the CNfix configuration), the leaf C/N ratio was set to the time-average PFT-average which was also inferred from the pd-CNdyn reference simulations.

10
The sensitivity of the simulated GPP to [CO 2  twice its present-day value along the 100 years. Both set-ups were repeated for the CNdyn and CNfix-time configurations, thus, resulting in a total of four idealized tests (1%CO2-CNdyn, 1%CO2-CNfix-time, 2xCO2-CNdyn and 2xCO2-CNfixtime simulations, see Table 2). For these sensitivity tests, climate drivers and nitrogen deposition files corresponding to the period of in-situ observations were used cyclic.

20
At global scale, in addition to the reference S1-CNdyn simulation, sensitivity tests were set-up to disentangle the main drivers of the simulated global increase in GPP (see Table 2). Based on the N 2 O Model Intercomparison Project protocol (Tian et al., 2018), additive scenarios were developed where GPP is only driven by climate change (S6-CNdyn); climate change and [CO 2 ] increase (S5-CNdyn); climate change, [CO 2 ] increase and land use change (S4-CNdyn); climate change, [CO 2 ] increase, land use change and nitrogen deposition evolution (S3-CNdyn); climate change, [CO 2 ] increase, land use 25 change, nitrogen deposition and nitrogen fertilizer evolution (S2-CNdyn) using the different forcing data presented in Sect.
2.3. A control simulation named S0-CNdyn was also performed, which extends the spin-up simulation (i.e., using the nitrogen input data, the [CO 2 ] value and land-use of the year 1860 and recycling the meteorological data of the period 1901-1920).

30
To test the impact on the GPP evolution of accounting for the carbon-nitrogen interactions, an additional simulation was run at global scale in which the gridded leaf C/N ratio was fixed to the time-average value for each PFT, inferred from the S0-CNdyn control simulation. Thus, the plant nitrogen status of each PFT within each grid-cell is the one of the pre-industrial era (i.e. for the year 1860). In this simulation, all drivers were varying over the period 1861-2016. However, the GPP being un-sensitive to variation of nitrogen deposition, nitrogen fertilization, and N from manure in the CNfix configuration, this simulation corresponds to a S4 scenario and is consequently named S4-CNfix.

Evaluation metrics
The mean seasonal cycle of the simulated GPP was evaluated against observations. This was done at the PFT level by 5 computing the mean seasonal cycle, averaged for all years and all sites belonging to a PFT. Alternatively, simulated daily GPP fluxes were evaluated by computing their root mean squared deviation (RMSD) to the observations. Further, the mean squared deviation (MSD) of the daily fluxes to the observation was decomposed and the contribution to the overall MSD from the mean bias, standard deviation and correlation was quantified based on the method of .
MSD can be written as: 10 where SB is the squared bias, SDSD is the squared difference between standard deviations and LCS the lack of correlation weighted by the standard deviations. SB, SDSD and LCS reflect respectively errors on bias, standard deviation and correlation. Finally, annual mean GPP flux was computed at each Fluxnet site in order to evaluate the model capacity at simulating site-to-site variations within each PFT. 15 The global S1-CNdyn simulation was also evaluated by comparing the spatial distribution of the annual mean GPP over the Lastly, the relative contributions of the different drivers to the present-day GPP value were assessed by computing the successive differences between additive scenarios. Thus, the contributions of climate change, [CO 2 ] increase, land use change, nitrogen deposition, nitrogen fertilization and nitrogen manure were provided respectively by the differences between S6 and S0, S5 and S6, S4 and S5, S3 and S4, S2 and S3, S1 and S2. Using this methodology, the sum of the 25 individual contributions is additive and equals the present-day GPP, thus ignoring possible non-linear interactions between drivers.

Evaluation of the standard configuration
The mean seasonal cycle of the GPP averaged per Plant Functional Type was reproduced for most PFTs, in the pd-CNdyn simulations (Fig. 1a). Over Tropical Evergreen Broadleaved Forests, the observed mean GPP values do not show a marked 5 seasonal cycle, while the simulated GPP are 3.0 gC m -2 day -1 higher in May-June compared to the rest of the year. For Temperate Evergreen Needleleaved Forests, the mean seasonal cycle simulated by ORCHIDEE matched the observations in terms of correlation (correlation coefficient of 0.99) and amplitude (model standard deviation of 2.4 gC m -2 day -1 to compare to 2.1 for observations), although the simulated GPP was overestimated by ~30% all the year. The model has the weakest performance for Temperate Evergreen Broadleaved Forest shown by a too pronounced seasonal cycle compared to that 10 observed (model standard deviation of 2.1 gC m -2 day -1 to compare to 0.7 gC m -2 day -1 for observations). In winter and early spring, GPP was overestimated by 2.0 gC m -2 day -1 , while later in June and July the decrease of the GPP was overestimated. The RMSE of the daily GPP flux averaged per PFT does not exceed 2.5 gC m -2 day -1 (Fig. 1b). The annual productivity of the PFTs being significantly different, the RMSE, when expressed as a percentage of the mean annual GPP (NRMSE), varies from 25% for Tropical Evergreen Broadleaved Forests to 80% for Boreal Evergreen Needleleaved Forests. RMSE was plotted against the length of the in-situ meteorological data (not shown) to check whether under-sampling of the climate 25 variability explained part of the bias. The data did not support such a relationship (correlation -0.218, CI 95% -0.448 -0.004). Figure 1c shows the respective contributions of SB (bias), SDSD (deviation) and LCS (correlation) on the total MSE per PFT. These relative contributions to the MSE differed depending of the PFT. Error on the mean bias was the largest contribution to MSE for Temperate and Boreal Evergreen Needleleaved Forests. At Tropical Evergreen Broadleaved Forests and C3 grassland sites, errors on the correlation contributed the most to the MSE while at Temperate Evergreen and 30 Deciduous Broadleaved Forest sites, the three sources of error were more equally distributed.
The simulated annual mean GPP per site was comparable to the one observed for most sites (Fig. 1d), with a RMSE averaged per PFT, which varied from 409 gC m -2 yr -1 (for Tropical Evergreen Broadleaved Forests) to 759 gC m -2 yr -1 (for Temperate Evergreen Broadleaved Forests). Both overestimation as well as underestimation were observed in all PFTs suggesting small systematic biases. Site-to-site variations of the annual mean GPP were relatively well reproduced for Tropical Evergreen Broadleaved Forests, Temperate Evergreen Needleleaved and Broadleaved Forests and C3 grassland 5 sites (Pearson's correlation coefficient of 1.0 but for only two sites, 0.63, 0.44 and 0.82, respectively) but not for Temperate Deciduous Broadleaved Forest and Boreal Evergreen Needleleaved Forest sites (Fig. 1d). For Temperate Deciduous Broadleaved Forest sites, the model produces significantly larger site-to-site GPP differences than the observations. In order to analyse how ORCHIDEE r4999 performs compared to r3977 (reference version without the nitrogen cycle), we 10 evaluated the GPP simulated by ORCHIDEE r3977 against Fluxnet observations (Fig. B1). The model/data agreement for r3977 was comparable to the one for r4999 but slightly better. In particular, NRMSE of the simulated daily GPP flux (Fig.   B1b) and the RMSE of the simulated annual mean GPP ( fig. B1d)  the PFT level for r3977 was due to a narrower range of NRMSE values at site level (whisker boxes are narrower), indicating that the NRMSE was not systematically lower at all sites but only at some specific ones.
Because the GPP flux is intimately linked to the transpiration flux through stomatal control, a site-level evaluation of the pd-CNdyn simulation has been performed against site-level observations of the latent heat (LE) flux (Fig. C1), an energy flux to 20 which transpiration contributes to, as does the soil evaporation. Overall, the model performed better at simulating LE variations than variations in GPP. This was particularly true when looking at the NRMSE of the simulated daily flux, which never exceeded 50% as a mean average score at the PFT level for LE, while it went to values up to 75% for GPP for some PFTs (BoENF and GRAc3).

Sensitivity to model configurations 25
The leaf C/N ratio in the pd-CNdyn simulation varied substantially over time and/or from one site to another (Fig. 2c). The observed variation was partly driven by the different nitrogen deposition load ( Fig. 2a and b) as well as by differences in the simulated mineralisation and plant Nitrogen uptake (not shown). V cmax being directly related to the leaf nitrogen content, the variations of the leaf C/N ratio induced seasonal variations of the V cmax on the order of 0 to 30%, depending of the sites.

30
The mean and median values of the MSE of the simulated daily GPP obtained from the range of sites within a vegetation class did not change substantially depending on whether nitrogen dynamics were accounted for or were fixed over time (pd-CNdyn, pd-CNfix-time or pd-CNfix-timePFT) (Fig. 3a). This finding holds for the error measures when decomposed into bias, standard deviation and correlation (Figs 3b-d). One exception to this model behaviour which is common across configurations is for some C3 grassland sites, where MSE and all its components were higher in the pd-CNfix-time and pd-CNfix-timePFT simulations compared to the pd-CNdyn simulation (Fig. 3). For Tropical and Temperate Evergreen Broadleaved Forest classes, the CNfix-time simulation exhibits narrower ranges for MSE or for some of its components (SB, SDSD, or LCS) compared to the pd-CNdyn simulation. This highlights the fact that, for some of the Tropical and Temperate 5 Evergreen Broadleaved Forest sites, constraining the leaf C/N ratio in time improved the fit to the observed GPP, while the same constraint deteriorated the fit at some other sites within that PFT. Results obtained with the pd-CNfix-timePFT simulation also lead to slightly narrower ranges of MSE values and of its bias (SB) subcomponent for Temperate Evergreen Needleleaved Forest, compared to the pd-CNfix-time simulation.

10
Simulated site-to-site variations of the annual mean GPP were sensitive to the model configuration regarding the leaf C/N ratio (Fig. 4). For all PFTs except the Boreal Evergreen Needleleaved Forest sites, the pd-CNdyn is the simulation, which best matched the observations in terms of RMSE. Nevertheless, the RMSE were of the same order of magnitude for the three model configurations for all PFTs except for C3 grassland sites for which the RMSE was significantly lower in the pd-CNdyn simulation (Fig. 4). 15

Sensitivity to [CO 2 ] increase
Increasing atmospheric [CO 2 ] by 1% per year lead to a continuous increase of the GPP in any of the two configurations; one with a dynamic leaf C/N ratios (1%CO2-CNdyn) and the other with leaf CN ratios fixed over time  and for any PFT class (Fig. 5). Note first that the large temporal cycle for the Tropical Evergreen Broadleaved forest is due to the recycling of two to four years of in situ meteorological forcing. As expected, the higher the [CO 2 ], the higher the GPP 20 was. Nevertheless, the GPP increase was less sensitive to the [CO 2 ] increase in the configuration with a dynamic C/N ratio (1%CO2-CNdyn), which reflects an induced Nitrogen-limitation of the photosynthesis. Notably, the GPP sensitivity to CO 2 increase in the 1%CO2-CNdyn simulation was particularly low for Boreal Evergreen Needleleaved Forest class. After 100 years of [CO 2 ] increase, the difference in GPP increase between the 1%CO2-CNfix-time and 1%CO2-CNdyn simulations reached 1.0 kgC m -2 yr -1 for Tropical Evergreen Broadleaved Forests 0.8 kgC m -2 yr -1 for Temperate Evergreen Needleleaved 25 Forest, 0.8 kgC m -2 yr -1 for Temperate Evergreen Broadleaved Forest, 0.6 kgC m -2 yr -1 for Temperate Deciduous Broadleaved Forest, 0.6 kgC m -2 yr -1 for Boreal Evergreen Needleleaved Forest and 1.0 kgC m -2 yr -1 for C3 natural grasslands. These differences in GPP increase corresponded to values of the order of 30-50% of the present-day annual mean GPP.
In the simulation where [CO 2 ] was doubled compared to the present-day level but N-limitation was not accounted for 30 (2xCO2-CNfix-time), GPP increased by an overall 90% (1.1 kgC m -2 yr -1 , Fig. 5), and at some sites even by as much as 150%. Accounting for a N-limitation (2xCO2-CNdyn simulation) reduced the overall GPP increase to ~50% compared to present-day value (0.6 kgC m -2 yr -1 , Fig. 5

). A decreasing trend in GPP increase for all PFTs -except Temperate Deciduous
Broadleaved Forest and C3 grassland sites -was apparent for the 2xCO2-CNdyn simulation, while such a trend was absent in the simulations without carbon-nitrogen interactions (2xCO2-CNfix-time). For instance, mean GPP increase at Temperate Evergreen Needleleaved Forest sites reached 0.8 kgC m -2 yr -1 over the 5 years consecutive to the doubling of the CO2 level but was only 0.5 kgC m -2 yr -1 after 100 years (Fig. 5). Furthermore, when comparing the simulations without and with Nlimitations, the year-to-year variability in GPP was amplified compared to the present-day variability for the configuration 5 without N-limitation. GPP increase induced by increase of atmospheric [CO 2 ] was achieved at a limited water cost or even resulted in saving water, at most sites. In the 1%CO2-CNfix-time simulations, after 100 years of [CO 2 ] increase, transpiration rate averaged per PFT was lower by few millimetres up to 50 mm per year, compared to its value in pd-CNfix-time simulation (Fig. 6, left  10 column). Averaged over all PFT, the mean transpiration rate decreased by ~25 mm yr -1 after 100 years. In the 1%CO2-CNdyn simulation, the decrease of the transpiration rate was even stronger (Fig. 6). Mean transpiration rate decrease averaged per PFT reached 75 mm to 110 mm per year, after 100 years. Averaged over all sites, the mean decrease equalled 100 mm per year, which corresponds to ~20% of its value in the pd-CNdyn simulation. Similar model behaviour with the 2xCO2 type simulations was exhibited (Fig. 6, right column). Compared to its value in the pd-CNfix-time simulation, the 15 mean transpiration rate averaged over all sites was stable in the 2xCO2-CNfix-time simulation. In the 2xCO2-CNdyn simulation, the mean transpiration rate was 50 mm per year lower than it was in the pd-CNdyn simulation (15%).
When expressed in terms of water use efficiency (WUE, unitless) -defined here as the ratio of GPP (gC m -2 yr -1 ) to transpiration (gH 2 O m -2 yr -1 ) -these model responses translated into an increase of the WUE for all the sensitivity tests (Fig.  20   7). Under the 1%CO2-CNfix-time simulation, after 100 years, WUE increased by 120% on average per PFT (Fig. 7, mean increase of 120%) compared to the pd-CNfix-time simulation. After 100 years, the mean WUE increase in the 1%CO2-CNdyn simulation is 17% lower than in the 1%CO2-CNfix-time simulation. Similarly, the mean WUE averaged per PFT in the 2xCO2-CNfix-time simulation is 70% to 90% higher than in the pd-CNfix-time simulation (mean increase of 80% over all sites, Fig. 7). In the 2xCO2-CNdyn simulation, the mean increase of the WUE is 10% lower than in the 2xCO2-CNfix-25 time simulation.

Global scale simulations
The present-day annual mean GPP under the S1-CNdyn simulation reached its maximum values in Central Africa (Fig. 8a).
Compared to values reported by MTE-GPP (Jung et al. 2010), simulated annual mean GPP in this region is overestimated by up to 2.0 kgC m -2 yr -1 (Fig. 8b). On the other hand, annual mean GPP values in the Amazon region or Indonesia were 30 underestimated by 1.5 kgC m -2 yr -1 . Overall, the Root Mean Square Deviation to the MTE-GPP values equals 0.7 kgC m -2 yr -and MTE-GPP data being less than 0.5 kgC m -2 yr -1 in most extra-tropical regions, the Root Mean Square Deviation is thus dominated by the mismatch between the simulations and observations over the tropics.
When summed by latitudinal band, the overlap between simulated annual mean GPP and MTE-GPP estimates increases ( Fig.   9) compared to the overall global mapping result. From 1982 to 2008, the simulated annual mean GPP above 25°N increased 5 from 36 PgC yr -1 to 42 PgC yr -1 (around 17% increase), with large year-to-year variations (up to 3 PgC yr -1 ) (Fig. 9a). Over the same period and the same domain, the MTE-GPP estimates varied between 37 and 40 PgC yr -1 (around 8% increase) thus showing a weaker positive trend compared to the ORCHIDEE simulation. In the tropical regions (between 25°S and 25°N), the MTE-GPP estimate slightly increased from 72 PgC yr -1 to 75 PgC yr -1 (around 4% increase) while according to the ORCHIDEE S1-CNdyn simulation, GPP between 1982 and 2008 varied from 62 PgC yr -1 to 72 PgC yr -1 (around 16% 10 increase) (Fig. 9b). In the tropical regions, the ORCHIDEE GPP strongly increases between 1992 and 2000 (from 64 to 72 PgC yr -1 ) and then fluctuates between 69 and 72 PgC yr -1 after 2000.
Due to a lower land area and a lower productivity per unit of land, the contribution of the southern lands (below 25°S) to the global GPP is rather small. It amounts 5 PgC yr -1 with ORCHIDEE, and 6 PgC yr -1 with MTE-GPP over the 1982-2008 15 period (Fig. 9c). Both estimates have no trend over the time period and the year-to-year variations are slightly larger with the ORCHIDEE model (standard deviation of 0.3 PgC yr -1 compared to 0.2 PgC yr -1 with MTE-GPP). Globally, annual mean GPP reached ~120 PgC yr -1 in 2016 (Fig. 9d). Fig. D1a to Fig. 8a, and Fig. 20 D1b to Fig. 8b) suggest that the bias in GPP originates from the bias in LAI rather than from more fundamental issues with the calculation of GPP. The model/data agreement for LAI when averaged per latitudinal band (Fig. E1) is comparable to the one for GPP (Fig. 9), with a good agreement for the Northern and Tropical lands and model underestimation in the Southern lands.

25
The agreement between the modelled and observed annual mean LAI and GPP summed over three latitudinal bands as well as at the global scale was higher for r4999 (ie S1-CNdyn simulation) compared to r3977 without the nitrogen cycle (see Fig.   9 and E1). r3977 systematically overestimated LAI and GPP for any region, except for the Southern lands where r3977 provided similar values than the GIMMS and MTE-GPP products, respectively. Compared to GIMMS and MTE-GPP products, gridded annual mean LAI and GPP values simulated by r3977 were overestimated in the Northern lands with 30 biases exceeding those found in r4999. On the opposite, biases of the r4999 were higher than those of r3977 in the tropical regions, in particular in Central Africa (see Fig. 8 and D1).

Discussion
We presented revision 4999 of the ORCHIDEE model that accounts for carbon-nitrogen interactions based on the initial 15 work of . In this model version, improvements were implemented with respect to water, energy and carbon budgets as well as to the nitrogen cycle and its coupling to the carbon cycle, compared to the original version of ORCHIDEE used in . Among these changes, a new carbon assimilation scheme was introduced  in which the maximum Rubisco activity-limited carboxylation rate is a direct function of the leaf nitrogen content ). Because these model developments primary impact on the GPP and because this flux 20 is the primary carbon input flow into the ecosystem, we focused our study on an in-depth evaluation of the simulated GPP flux, from site to global scale, and from daily to mean annual time scale.
We showed that the version of ORCHIDEE (r4999) that accounts for carbon-nitrogen interactions can simulate reasonably well the mean seasonal cycle, daily mean variations, and annual mean GPP for most PFTs. The overlap between the 25 ORCHIDEE simulations and the eddy-covariance based observations is similar to other models Slevin et al., 2015). For instance, based on an evaluation over 32 Fluxnet sites,  reported mean RMSE values of the daily GPP flux simulated by three GTEMs (ie. ORCHIDEE, CTESSEL and ISBA-A-gs) which range between 1.9 and 4.4 gC m -2 d -1 depending of the model and PFT considered (excluding cropland sites for which all models performed the least and which are not included in our study). In comparison, mean RMSE of the simulated daily GPP flux in our study 30 for any PFT never exceeded 2.5 gC m -2 d -1 .
The exception to the general agreement is for the Temperate Evergreen Broadleaved Forest sites for which the ORCHIDEE r4999 failed to reproduce the observed mean seasonal cycle. Indeed, the weak performances are for the five Mediterranean sites, but not for the two Australian sites. The Mediterranean sites receive high N-deposition loads leading to low leaf C/N ratios. Over these sites, at the beginning of the growing season, the higher maximum Rubisco activity-limited carboxylation rate values (Vc max ) induced by high leaf nitrogen content is the main reason explaining the GPP overestimation. The 5 overestimation of the GPP at the early stage of the growing season induced higher rates of transpiration. These higher rates of transpiration partly explain the underestimation of the GPP during the summer, due to a too strong depletion of the soil water. Because the too low GPP during the summer demands less nitrogen to support biomass growth, this results in higher nitrogen available later in the season. This mechanism tends to maintain or even amplify the mismatch at the beginning of the growing season. The two Australian sites receive low N-deposition loads in comparison to the three Mediterranean sites. 10 Consequently, the Australian sites have high leaf C/N ratios (see Fig. 2c, for TeEBF where Australian sites correspond to the blue and green lines), which overcomes the overestimate of the GPP at the beginning of the growing season and the subsequent issues. Nevertheless, GPP appears significantly biased -against MTE-GPP -over the tropical regions, with positive biases in Central Africa and negative ones in Amazonia. One may note that these biases are not so contrasted in the original ORCHIDEE version without the nitrogen cycle (r3977, see Figure 8c). Further analyses showed that the GPP biases over tropical forest regions are driven by different leaf C/N ratios across regions. However, it remains unclear what are the 25 primary drivers of the spatial variation of the leaf C/N ratio and consequently of GPP. One of the drivers is likely to be NOx deposition, which is lower in Amazonia compared to Central Africa (not shown). There is no such relationship between GPP and BNF rate, nor between GPP and NHx deposition rate, in tropical regions. The drivers and/or processes that are responsible for turning large scale differences in NOx deposition into fine-scale differences in GPP are yet to be identified.

30
Apart from the MTE-GPP product, there is few spatially explicit data available that is suitable to evaluate global simulations of GPP. Due to the fact that two CO 2 fluxes occur at land but have an opposite direction (ie GPP and Total Ecosystem Respiration) data of atmospheric [CO 2 ] cannot be used for this purpose. As an alternative,  assessed the GPP growth rate over the 20 th century based on an independent estimate by using atmospheric carbonyl sulfide (COS) records. Based on a procedure that optimizes the COS sources and sinks to match the observed atmospheric [COS] dynamic, they deduce that the GPP increased by 31%±5% over the 20 th century. This corresponds to a larger GPP growth rate than that estimated by most Global Terrestrial Ecosystem Models. Although higher than the 31%±5% based on the COS constraint, the ORCHIDEE based of 39% is in good agreement with the COS constraint compared to most of other ecosystem models reported in the study of . This increases our confidence in the ability of ORCHIDEE r4999 to account 5 for the interactions between carbon and nitrogen.
Constraining the leaf C/N ratio over time and per PFT (using the values from a simulation with dynamic C/N ratios) was found not to impact model performance at site-level, for most PFTs. Note, however, that the leaf C/N ratios in the pd-CNfixtime and pd-CNfix-timePFT simulations were set using the present-day values from the pd-CNdyn simulation. At some sites, 10 the mean standard error of the simulated daily mean GPP variation is reduced, but an increase in mean standard error was observed at other sites. At first, the similarity in performance for configurations with and without dynamic leaf C/N ratio may appear as a failure in term of expected impacts with the inclusion of the nitrogen cycle in the model. Accounting for dynamic leaf C/N ratio and for site-specific information such as nitrogen deposition could be expected, at least from a theoretical point of view, to enhance the predictive power of ORCHIDEE. However, at the same time, dynamic modelling of 15 the leaf C/N ratio adds complexity and interactions, which come with additional sources of uncertainties. Our results indicate that for some sites the increase in predictive power dominates whereas for other sites the additional uncertainties dominate the total model error.
One exception were the C3 grassland sites for which the mean model performance was lower when carbon-nitrogen 20 interactions were not accounted for (pd-CNfix-time and pd-CNfix-timePFT simulations), leading to a higher mean MSD for this PFT compared to the pd-CNdyn simulation. Constraining or not the leaf C/N ratio (pd-CNfix-time vs. pd-CNdyn configuration) alters GPP via changes on the Vc max value. Moreover, for some "extreme" conditions in the CNdyn configuration, the carbon allocation to the different reservoirs may be reduced due to insufficient nitrogen in the labile pool, irrespective of the value of the C/N ratio of the standing leaf biomass (see end of Sect. 2.1.3 where this model behaviour is 25 detailed). These extreme conditions and their consequences on the biomass allocation cannot be captured with the CN-fix configuration and it appears that at many C3 grassland sites, these extreme conditions are sufficiently frequent to produce different model responses and performances between the pd-CNfix-time and pd-CNfix-timePFT simulations and the pd-CNdyn simulation.

30
While behaviour and performance are similar across ORCHIDEE configurations for present-day conditions for most PFTs, they differ substantially in terms of the response of GPP to enhance atmospheric [CO 2 ]. Site-mean increase in GPP after 100 years under the 1%CO2 experiment is half when carbon-nitrogen interactions are accounted for compared to ignoring the carbon-nitrogen interactions (CNdyn, compared to CNfix-time). Interestingly, this effect-size is of the same order of magnitude as the simulated reduction in the increase in NPP (65%) under the Representative Concentration Pathway 8.5  estimated by coupled carbon-climate model projections and C/N stoichiometric models. The reduction in the increase in NPP was caused by the interaction between enhanced atmospheric [CO 2 ] and nitrogen limitations trends .

5
Simulated site-average GPP increased by 50% when [CO 2 ] was doubled (~+370 ppm, 2xCO2 experiment) and carbonnitrogen interactions were accounted for. A similar response has been reported for some experimental FACE sites  where GPP increased by 30% for an increase between 150 and 200 ppm in atmospheric [CO 2 ].
Also, the simulated decrease of the growth rate in GPP over time at some sites is in line with the Progressive Nitrogen Limitation as postulated by . At global scale, we showed that the present-day simulated GPP is reduced by 10 20% when the additional nitrogen limitation associated to the period 1860-2016 is accounted for (S4-CNdyn) compared to ignoring this nitrogen limitation (S4-CNfix). We showed that the simulated GPP trajectories from 1861 to 2016 under the S1-CNdyn and the S4-CNfix (where the leaf C/N ratios are fixed to their mean-time pre-industrial values) are similar, meaning that the global historical increase of nitrogen inputs (through deposition and fertilisation) was of the same order than the nitrogen demand to fulfil the GPP increase due to the [CO 2 ]-fertilisation effect. This general behaviour, however, 15 hides contrasted responses between classes of PFTs. Over forested lands, the potential GPP increase primary driven by the [CO 2 ]-fertilisation effect has been significantly limited due to a shortage in plant available nitrogen (compare S4-CNfix and S3-CNdyn curves of fig. 10b). On the opposite, over grasslands and croplands, the supply of nitrogen through fertilisation has fostered the GPP more than it would have done solely from the [CO 2 ]-fertilisation effect without N-limitation (compare S4-CNfix and S1-CNdyn curves of fig. 10c). 20 We showed that accounting for carbon-nitrogen interactions or not, when atmospheric [CO2] is enriched does not only impact on the GPP response but has also important consequences on the transpiration rate. We showed that on average, when atmospheric [CO 2 ] is doubled (~ 700-750 ppm), the water use efficiency -calculated as the ratio between GPP and transpiration -increases by 80% under the CNfix-time configuration, but only by 67% under the CNdyn configuration. The 25 difference in water use efficiency change is likely the outcome of partly the nitrogen cycle through its impact on Vc max and thus on stomatal conductance. Although based on a slightly different indicator, this is well in agreement with the observed increase (+68%) of the instantaneous transpiration efficiency (computed as the assimilation divided by the stomatal conductance, ITE) reported in the meta-analysis of  based on a set of C3 ecosystem sites where atmospheric [CO2] was enriched to reach atmospheric concentrations between 475 and 600 ppm. Our model experiments 30 illustrate the interplay between carbon, nitrogen and water cycles and highlight the need for careful analysis when modelling ecosystem productivity under future climate changes and its possible reduction due to water limitation .

Concluding remarks
We conclude that if carbon-nitrogen interactions are accounted for, the functional responses of ORCHIDEE r4999 better agrees with current understanding of photosynthesis than when the carbon-nitrogen interactions are not accounted for. From this point of view our factorial experiment and sensitivity analysis confirm what has been shown by several other Global Terrestrial Models, i.e., that carbon-nitrogen interactions are essential in understanding global terrestrial ecosystem 5 productivity (Goll et al., 2017;Sokolov et al., 2008;Thornton et al., 2009;Wang and Houlton, 2009;Wania et al., 2012; and especially its dynamic under atmospheric [CO 2 ] increase and climate change. Further simulations, making use of this new ORCHIDEE version and driven by different socio-economic scenarios over the 21 st century (i.e. Representative Concentration Pathways) will be performed in order to better quantify the future GPP response to the combined evolution of [CO 2 ] and nitrogen land supply. Additionally, climate simulations with the IPSL 10 Earth System Model (as part of the CMIP6 inter-comparison project) including this new ORCHIDEE version are planed to assess the impact of nitrogen limitation on the fate of the net land carbon sink and on climate projections.

Code availability
The source code is freely available online via the following address: http://forge.ipsl.jussieu.fr/orchidee/wiki/GroupActivities/CodeAvalaibilityPublication/ORCHIDEE_gmd-2018-261. A DOI 15 has been requested for this page which provides guidelines for downloading.

Author contribution
NV and PM developed the model code with contributions from SL, SZ, BG and JG. NV performed the simulations. NV and VB prepared the figures. NV prepared the manuscript with contributions from all co-authors.