Global evaluation of the nutrient-enabled version of the land surface model ORCHIDEE-CNP v1.2 (r5986)

The availability of phosphorus (P) and nitrogen (N) constrains the ability of ecosystems to use resources such as light, water and carbon. In turn, nutrients impact the distribution of productivity, ecosystem carbon turnovers and their net exchange of CO₂ with the atmosphere in response to variation of environmental conditions in both space and time. In this study, we evaluated the performance of the global version of the land surface model ORCHIDEE-CNP (v1.2), which explicitly simulates N and P biogeochemistry in terrestrial ecosystems coupled with carbon, water and energy transfers. We used data from remote sensing, ground-based measurement networks and ecological databases. Components of the N and P cycle at different levels of aggregation (from local to global) are in good agreement with data-driven estimates. When integrated for the period 1850 to 2017 forced with variable climate, rising CO₂ and land use change, we show that ORCHIDEE-CNP underestimates the land carbon sink in the Northern Hemisphere (NH) during recent decades despite an a priori realistic gross primary productivity (GPP) response to rising CO₂. This result suggests either that processes other than CO₂ fertilization, which are omitted in ORCHIDEE-CNP such as changes in biomass turnover, are predominant drivers of the northern land sink and/or that the model parameterizations produce emerging nutrient limitations on biomass growth that are too strict in northern areas. In line with the latter, we identified biases in the simulated large-scale patterns of leaf and soil stoichiometry as well as plant P use efficiency, pointing towards P limitations that are too severe towards the poles. Based on our analysis of ecosystem resource use efficiencies and nutrient cycling, we propose ways to address the model biases by giving priority to better representing processes of soil organic P mineralization and soil inorganic P transformation, followed by refining the biomass production efficiency under increasing atmospheric CO₂, phenology dynamics and canopy light absorption.

1 Introduction

Nitrogen (N) and phosphorus (P) are key macronutrients that control metabolic processes and plant growth and constrain ecosystem-level productivity (Elser et al., 2007; Norby et al., 2010; Cleveland et al., 2013). The amount and stability of soil carbon (C) stock are also affected by N and P through their regulating role in the mineralization of litter and soil organic matter (Gärdenäs et al., 2011; Melillo et al., 2011). The availability of N and P is likely to limit future carbon storage under climate change and rising atmospheric CO₂. Empirical stoichiometry observations were applied in the posteriori estimates of future carbon storage from land surface models (LSMs) lacking an explicit simulation of N and P biogeochemistry, which consistently led to an overestimation of future carbon storage in LSMs (Hungate et al., 2003; Wang and Houlton, 2009; Zaehle et al., 2015; Wieder et al., 2015). Nevertheless, this approach has large uncertainties (Penuelas et al., 2013; Sun et al., 2017) and relies on unproven assumptions (Brovkin and Goll, 2015).

An alternative is to directly represent the complex interactions between N, P and carbon in an LSM. Several LSMs incorporated different parameterizations of N interactions (e.g., Thornton et al., 2007; Zaehle et al., 2014), but very few global models have included P interactions. The few models accounting for P limitation in plant growth showed that P availability limits primary productivity and carbon stocks on highly weathered soils in the tropics (Wang et al., 2010; Yang et al., 2014), and one study also suggested that P limitations could also occur in the Northern Hemisphere in the near future (Goll et al., 2012). Model representations of P interactions are highly uncertain since the critical processes are poorly constrained by current observational data. In particular, the desorption of P from soil mineral surfaces and the enhancement of P availability for plants by phosphatase enzymes secreted by plant roots and microbes were identified to be critical but poorly constrained (Fleischer et al., 2019).

Previous studies (Wang et al., 2010; Goll et al., 2012; Yang et al., 2014; Thum et al., 2019) have suggested that the inclusion of the phosphorus cycle improves model performances with regard to reproducing observed C fluxes. But adding new and uncertain P-related processes does not grant an automatic improvement in an LSM in general. First, more (nutrient-related) equations with more uncertain parameters can result in less robust predictions. Second, models ignoring nutrients were often calibrated on available carbon data so that a new model with nutrients inevitably needs a parameter recalibration to reach similar performances as the same model without nutrients. Third, for evaluating a large-scale model resolving both nutrient and carbon biogeochemistry, one needs specific nutrient-related datasets, which are more scarce than classical biomass, productivity and soil carbon data used for benchmarking carbon-only models.

The evaluation for N and P, together with carbon cycling in global LSMs, remains very limited (Wang et al., 2010; Goll et al., 2012), but recent advances in ground-based measurements, ecological datasets and process understanding have made a better evaluation of C, N and P models feasible. The available nutrient datasets have allowed for meta-analyses of site-level nutrient fertilization experiments (e.g., Yuan and Chen, 2015; Wright, 2019), data-driven assimilation schemes to constrain nutrient budgets (Wang et al., 2018), new knowledge about the critical P processes of sorption (Helfenstein et al., 2018, 2020) and phosphatase-mediated mineralization (Sun et al., 2020), global datasets of leaf nutrient content (Butler et al., 2017), and empirical constraints on the CO₂ fertilization effect on land carbon storage (Terrer et al., 2019; Liu et al., 2019). In addition to direct comparison with nutrient datasets, it is also possible to diagnose emerging model responses in terms of ecosystem resource use efficiencies (RUE) and confront them with observations for identifying how ecosystems adjust and optimize nutrient, water, light and carbon resource availabilities (Fernández-Martínez et al., 2014; Hodapp et al., 2019). In particular, modeled N and P use efficiencies can be compared to observation-based estimates at ecosystem scale (Gill and Finzi, 2016) and at biome scale (Wang et al., 2018).

Here we evaluate the global cycles of C, N and P in the nutrient-enabled version of the LSM ORCHIDEE, ORCHIDEE-CNP (v1.2). The model has been previously evaluated for tropical sites (Goll et al., 2017a, 2018) and for coarse-scale global carbon fluxes and stocks using the International Land Model Benchmarking system iLAMB by, e.g., Friedlingstein et al. (2019). The results from this evaluation showed a slightly worse performance for ORCHIDEE-CNP (v.1.2) than the carbon-only version of ORCHIDEE, which has been extensively calibrated (Friedlingstein et al., 2019). In this study, we perform a detailed evaluation of ORCHIDEE-CNP focusing on four ecosystem characteristics that were found to be critical for the response of land C cycling to increasing CO₂ and climate change: (1) vegetation resource use efficiencies, (2) the response of plant productivity to increasing CO₂, (3) ecosystem N and P turnover and openness, and (4) large-scale patterns of ecosystem stoichiometries. Points (1) and (2) control the response of vegetation carbon storage operating on timescales of years to decades, while points (3) and (4) control the carbon storage potential on an ecosystem level, which determines the response on much longer timescales. Further, the implications of including nutrient cycles on the simulated land C cycling are discussed.

Figure 1 Schematic of C, N and P cycles considered in ORCHIDEE-CNP.

2 Modeling

2.1 Model description

ORCHIDEE-CNP simulates the exchange of greenhouse gases (i.e., carbon dioxide, nitrous oxide), water and energy at the land surface and features a detailed representation of the root uptake of dissolved N and P, the allocation of N and P among tissues, and the N and P turnover in litter and soil organic matter (Goll et al., 2017a, 2018) (Fig. 1). In this study, we present the first global application of the model and an evaluation against global carbon and nutrient datasets. ORCHIDEE-CNP simulates the cycles of C, N and P, which are described in detail elsewhere (Krinner et al., 2005; Zaehle and Friend, 2010; Goll et al., 2014, 2017a, 2018). We give a brief overview here. P enters the ecosystem through release from minerals into the soil solution, whereas N is biologically fixed from an ample reservoir of dinitrogen. Dissolved nutrients are either taken up by vegetation, converted into soil organic matter or absorbed onto soil particles. Losses occur as leaching of dissolved nutrients, gaseous soil N emissions or occlusion of P in secondary minerals. When nutrients are taken up by vegetation they are either stored internally or used to build new plant tissue driven by the availability of C, N and P in vegetation. The nutrient concentration of plant tissue varies within a prescribed range depending on the relative availability of C, N and P. Before plant tissue is shed, depending on the tissue a fixed fraction of the nutrients is recycled. The nutrients contained in dead plant tissue and organic matter are mineralized and released back into the soil solution. The model version applied in this study is based on Goll et al. (2017a, 2018) and referred to as ORCHIDEE-CNP v1.2. Major modifications compared to v1.1 are described as follows (details can be found in Sect. S1 in the Supplement).

The original formulation of photosynthetic capacity in ORCHIDEE-CNP v1.1 assumed leaf N to be the sole regulator of leaf photosynthetic characteristics (Kattge et al., 2009). Here, we applied a new empirical function that relates photosynthetic capacity to both the leaf N and P concentration based on data from 451 species from 83 different plant families (D. S. Ellsworth et al., unpublished data). Leaf C : N : P ratios that were prescribed in ORCHIDEE-CNP v1.1 a priori in a narrow range specific to the plant functional type (PFT) are now given a larger range common to all PFTs (Table S1), allowing for the prediction of variation in leaf stoichiometry across climate and soil gradients independently of the prescribed vegetation (PFT) map.

In ORCHIDEE-CNP v1.1, an empirical function, f(T_soil), was used to reduce biochemical mineralization and plant nutrient uptake at low soil temperature (Eq. 5 in Goll et al., 2017a), which was adopted from the N-enabled version of ORCHIDEE (Zaehle and Friend, 2010) to avoid an unrealistic accumulation of N within plants when temperatures are low. We found that this function was not needed when P uptake is accounted for and was thus removed. It should be noted that this temperature dependence is different from the one that describes the temperature dependence of soil organic matter (SOM) and litter decomposition. For grasslands and croplands, we implemented root dormancy, which is triggered by drought or low temperatures. During dormancy, root maintenance respiration is reduced by 90 % following Shane et al. (2009), but root acquisition of soil nutrients continues as long as root biomass exists (Malyshev and Henry, 2012). It should be noted that total root loss can occur for extremely long droughts or cold periods when maintenance respiration depletes root carbon.

Several parameters were recalibrated, i.e., the coefficient relating maintenance respiration to biomass and the leaf to sapwood ratio, or corrected in the case of the turnover of sapwood for tropical evergreen broadleaf forest (TREBF) and tropical rain-green broadleaf forest (TRDBF) to achieve more realistic wood growth rates for those forests (not shown). We also adjusted the recycling efficiency of nutrients from the root ($f_{\text{trans,\hspace{0.17em}root}}^{\text{N}}$, $f_{\text{trans,\hspace{0.17em}root}}^{\text{P}}$) and leaf ($f_{\text{trans,\hspace{0.17em}leaf}}^{\text{N}}$, $f_{\text{trans,\hspace{0.17em}leaf}}^{\text{P}}$) according to data compilations from Freschet et al. (2010) and Vergutz et al. (2012). The new values of these parameters and their sources are given in Supplement (Sect. S1).

2.2 Simulation setup

We performed a global simulation at 2^∘ × 2^∘ spatial resolution for the historical period (1700–2017) by adapting the TRENDY version 6 protocol (Sitch et al., 2015; Le Quéré et al., 2018). The simulation was performed using historical climate forcing, land cover changes and management (i.e., mineral fertilizer application, crop harvest; see Sect. 3.1.6), and atmospheric CO₂ concentrations (S3-type simulation). Prior to the historical simulation, we performed a model spin-up to equilibrate the C, N and P pools and fluxes (Sect. S1A in the Supplement) by forcing the model with cycled climate forcing of 1901–1920 and the land cover map and land management corresponding to the year 1700. To disentangle the effect of introducing nutrient cycles into ORCHIDEE, we performed the same simulation with ORCHIDEE (revision 5375), which has no nutrient cycles and a comparable parameterization for other processes. ORCHIDEE was run at a higher spatial resolution (0.5^∘ × 0.5^∘) than ORCHIDEE-CNP. Prior to the analysis, the data from ORCHIDEE were remapped to the resolution of ORCHIDEE-CNP.

2.2.1 Meteorological data

The model was forced by CRU-JRA-55 meteorological data provided at a spatial resolution of 0.5^∘ × 0.5^∘ and upscaled to a resolution of 2^∘ × 2^∘. These data comprise global 6-hourly climate forcing data providing observation-based temperature, precipitation and incoming surface radiation. They are derived from Climatic Research Unit (CRU) TS3.1 monthly data (Harris et al., 2014) and the Japanese 55-year Reanalysis (JRA-55) data (Kobayashi et al., 2015), covering the period 1901 to 2017. This climate dataset was provided by the TRENDY v6 model intercomparison project (Le Quéré et al., 2018).

2.2.2 Land cover

The historic land cover change maps were based on the European Space Agency Climate Change Initiative (ESA-CCI) land cover data (Bontemps et al., 2013). To be used by global vegetation models like ORCHIDEE-CNP, ESA-CCI land cover data were aggregated to 2^∘ × 2^∘ and grouped into PFTs using the reclassification method from Poulter et al. (2011, 2015). The fraction of cropland and pasture in the PFT map was further constrained by the cropland area and the sum of pasture and rangeland area for the year 2010, respectively, in the History Database of the Global Environment land use dataset (HYDE 3.2; Klein Goldewijk et al., 2017a, b), which were also aggregated to 2^∘ × 2^∘. The above processes produced a reference ESA-CCI-based PFT map for the year 2010. The land use changes derived from Land-Use Harmonization (LUH) v2 (http://luh.umd.edu/data.shtml, last access: 3 September 2018; an update release of Hurtt et al., 2011) were aggregated to 2^∘ × 2^∘ and then applied to this reference PFT map to constrain the land cover changes of forest, grassland, pasture and rangeland, and cropland during the period 1700–2017 using the backward natural land cover reconstruction method of Peng et al. (2017). As a result, a set of historic PFT maps suitable for global vegetation models was established, distinguishing global land cover changes for the period of 1700–2017 at 2^∘ × 2^∘ resolution.

2.2.3 Soil and lithology datasets

ORCHIDEE-CNP v1.2 is forced by information on soil texture, pH, bulk density and soil types (Goll et al., 2017a). We used a global gridded map of three soil texture classes from Zobler (1986) to derive soil-texture-specific parameters for soil water capacity, hydraulic conductivity and thermal conductivity. We used global gridded data on bulk density from the Harmonized World Soil Database (HWSD; FAO/IIASA/ISRIC/ISSCAS/JRC, 2012) and soil pH from the International Geosphere–Biosphere Programme Data Information System Soil Data (Global Soil Data Task Group, 2000). Soil pH forcing maps are needed to simulate the dynamics of NH₃ and NH${}_{\mathrm{4}}^{+}$ in soil in ORCHIDEE (Zaehle and Friend, 2010). We used a global gridded map with the dominant soil orders (following the USDA Soil Taxonomy) at 1^∘ × 1^∘ resolution to derive soil-order-specific soil phosphorus sorption parameters (Goll et al., 2017a).

The release of P from chemical weathering of rocks is computed dynamically following Goll et al. (2017a) and depends on the lithology types and soil shielding (discontinuation of the active soil zone from the bedrock) (Hartmann et al., 2014). We used the global lithological map (GLiM) of Hartmann and Moosdorf (2012) upscaled to 1^∘ × 1^∘ resolution, which accounts for the lithology fractional coverage of 16 classes on a sub-grid scale. We also used a spatially explicit map of soil shielding at 1^∘ × 1^∘ resolution (Hartmann et al., 2014).

2.2.4 Atmospheric nitrogen and phosphorus deposition

Global gridded monthly atmospheric N and P deposition during 1860–2017 was derived from a reconstruction based on the global aerosol chemistry–climate model LMDZ-INCA (Wang et al., 2017). LMDZ-INCA was driven by emission data, which included sea salt and dust for P, primary biogenic aerosol particles for P, oceanic emissions for N (NH₃), vegetation emissions for N (NO), agricultural activities (including fertilizer use and livestock) for N, and fuel combustion for both N (NO_y and NH_x) and P. Reconstructions for the years 1850, 1960, 1970, 1980 and 1990, as well as each year from 1997 to 2013, were linearly interpolated to derive a time series for 1850–2013. For the period before 1850, we assumed N and P deposition rates of the year 1850. For the period after 2013, we assumed rates of the year 2013. In ORCHIDEE-CNP, atmospheric N and P deposition is added to the respective soil mineral N and P pools without considering interception by the canopy.

2.2.5 Nutrient management

For croplands, we used yearly gridded mineral N and P fertilizer application data from Lu and Tian (2017) available for the period 1960 to 2017. This dataset is based on national-level data on crop-specific fertilizer application amounts from the International Fertilizer Industry Association (IFA) and the FAO. N and P mineral fertilization between 1900 and 1959 was linearly extrapolated assuming that fertilizer applications for 1900 are zero and that there were no N and P fertilizers applied before 1900. For pasture, we used global gridded datasets of mineral N fertilizer application rates from Lu and Tian (2017), developed by combining country-level statistics (FAO) and land use datasets (HYDE 3.2) (Xu et al., 2019). For both cropland and pasture, N and P in mineral fertilizer were assumed to go directly into soil mineral pools, where all mineral N fertilizer was assumed to be in the form of ammonium nitrate, with half of N as ammonium (NH${}_{\mathrm{4}}^{+}$) and half as nitrate (NO${}_{\mathrm{3}}^{-}$).

Manure applications are also included as a model forcing, given their significant input contribution to agricultural soils. For cropland, we used gridded annual manure N application data for the period 1860–2014 from Zhang et al. (2017) compiled and downscaled based on country-specific annual livestock population data from FAOSTAT. For the period before 1860, we assumed N and P deposition rates of the year 1860. For pasture, we used global gridded datasets of N manure application rates from Lu and Tian (2017). The application of manure P in cropland and pasture was derived from manure N assuming a manure P : N ratio of 0.2. This ratio is a weighted value by the amount of manure N applied to soil and derived from ruminants (14.4 Tg N yr⁻¹) and monogastric animals (10.1 Tg N yr⁻¹) from FAOSTAT for the year 2000 with P : N ratios of 0.165 in ruminant manure (mean of 0.15–0.18 from Lun et al., 2018) and 0.26 in monogastric manure (mean of 0.24–0.28 from Lun et al., 2018). For manure applied to cropland and pasture, we assumed a typical slurry application with 90 % of the N in the liquid part of the slurry (like urine) going into the soil NH4⁺ pool. For the solid part of the slurry, we assumed it goes into a litter pool with a C : N ratio of 10 : 1 following Soussana and Lemaire (2014).

Mineral and manure N and P fertilizers in cropland were applied at day of year (DOY) 120 for the Northern Hemisphere (30–90^∘ N), DOY 180 for tropical regions (30^∘ N–30^∘ S) and DOY 240 for the Southern Hemisphere (30–90^∘ S).

Table 1 Main information on datasets used for global evaluation of ORCHIDEE-CNP.

3 Evaluation

We evaluated the performance of ORCHIDEE-CNP v1.2 based on four major aspects (Fig. 1). Firstly, we evaluated the global C, N and P flows and storages. In the absence of robust spatially resolved estimates of N and P fluxes, we used the data-driven reconstruction of steady-state C, N and P fluxes on the biome level from the Global Observation-based Land-ecosystems Utilization Model of Carbon, Nitrogen and Phosphorus (GOLUM-CNP) v1.0 (Wang et al., 2018) (Table 1). Secondly, we evaluated plant resource use efficiencies (RUEs) of light, water, C, N and P on global and biome scales. RUEs reflect how ecosystems adjust and adapt to the availability of nutrient, water, light, and carbon resources (Fernández-Martínez et al., 2014; Hodapp et al., 2019). For this, we used estimates from site measurements and observation-based gridded datasets. Thirdly, we evaluated the response of gross primary productivity (GPP) to elevated CO₂ to assess the response of plant productivity to changing resource availability (i.e., CO₂) and historical perturbation C fluxes. For this, we used observation-based estimates (Ehlers et al., 2015; Campbell et al., 2017). Fourthly, we evaluated large-scale patterns of vegetation and soil N : P ratios as well as the N and P openness and turnover rates on the ecosystem level to assess spatial variation in nutrient limitation and the underlying drivers. For this, we used estimates from GOLUM-CNP, site measurements and observation-based gridded datasets (Kerkhoff et al., 2005; McGroddy et al., 2004; Reich and Oleksyn, 2004; Tipping et al., 2016; Butler et al., 2017; Wang et al., 2018). Finally, we showed the implications of ORCHIDEE-CNP for C cycling by evaluating the spatiotemporal patterns of terrestrial C fluxes and pools of the two versions of ORCHIDEE. For this, we used observation-based products of GPP and atmospheric inversions of the net land–atmosphere CO₂ flux excluding fossil fuel emissions (Table 1). Each dataset is summarized in Table 1 and described in detail in the Supplement. All the gridded datasets with high spatial resolutions (Table 1) were resampled to the 2^∘ × 2^∘ resolution of the model output using area-weighted mean methods.

3.1 Ecosystem productivity

Different data-driven maps of NPP and GPP based on remote sensing and climate datasets were used (Table 1), thereby accounting for the uncertainty of each product and for the uncertainty from the spread between different products. Uncertainties of each NPP and GPP product were derived according to the original publications. We used a 20 % uncertainty of gridded GPP from the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Breathing Earth System Simulator (BESS) (Sect. S1C in the Supplement; Turner et al., 2006; Jiang and Ryu, 2016) at a 2^∘ scale. This is a coarse extrapolation of uncertainty reported at the grid cell scale, since none of these products reported spatial error covariance information, allowing us to upscale this uncertainty at 2^∘ resolution. Further, for some products, uncertainty was defined as the bias against local measurements (Turner et al., 2006) and for others by using different climate input fields (Table 1). For multi-tree ensemble (MTE) GPP (Table 1), we used the spread (1σ standard deviation) from an ensemble of 25 members produced by different machine-learning methods (Jung et al., 2009). For MODIS NPP (Table 1), we used a 19 % uncertainty as assessed by Turner et al. (2006). For BETHY NPP we do not have an uncertainty (Tum et al., 2016). For Global Inventory Modeling and Mapping Studies (GIMMS) NPP (Table 1), we used the variance of three sets of products (Table 1) based on different climate datasets (Smith et al., 2016).

Two statistical indices were used to summarize the performance of ORCHIDEE and ORCHIDEE-CNP with respect to the interannual and seasonal variability of GPP and the interannual variability of net biome productivity (NBP) (Sect. 4.6): the coefficient of determination (R²) and relative mean square deviation (rMSE). The rMSE is defined as

$$\begin{array}{}\text{(1)}& \text{rMSE}={\displaystyle \frac{\sum _{j=\mathrm{1}}^n{\left(X_{\text{model},j}-X_{\text{ref},j}\right)}^{\mathrm{2}}}{\sum _{j=\mathrm{1}}^n{\left(X_{\text{ref},j}-{\widehat{X}}_{\text{ref},j}\right)}^{\mathrm{2}}}}.\end{array}$$

X_model and X_ref are values from models (i.e., ORCHIDEE and ORCHIDEE-CNP) and referenced datasets (i.e., MTE and BESS; Sect. S1C), respectively, and ${\widehat{X}}_{\text{ref},j}$ is the mean value across all years (for interannual variability evaluation) or all months (for seasonality evaluation).

3.2 Resource use efficiencies

The definition of resource use efficiencies is explained in Sect. 4.2. Observation-based light use efficiency (LUE) was calculated using MTE GPP, downward shortwave radiation from CRUJRA, and the fraction of absorbed photosynthetically active radiation (fAPAR) from the Global SeaWiFS Level-3 data (Gobron et al., 2006a, b). Uncertainty was derived from 25 ensemble members of MTE GPP. Observation-based water use efficiency (WUE) was calculated as the ratio between MTE GPP and MTE ET (Table 1); its uncertainties were calculated using a Monte Carlo resampling procedure in which 25 different members of GPP and ET were randomly selected. Observation-based carbon use efficiency (CUE) was calculated from the ratio of MODIS NPP to MODIS GPP. It should be noted that MODIS NPP is based on a calibrated version of the BIOME-BGC model (Turner et al., 2006) so that CUE is not strictly an observation-based quantity. CUE uncertainties were calculated using a Monte Carlo method given 20 % and 19 % uncertainty for MODIS GPP and NPP products at 2^∘ resolution, respectively.

Figure 2 Flowchart of mean flows and storages per area of C, N and P (g C g⁻¹, N g⁻¹, P m⁻² yr⁻¹) in natural biomes for GOLUM-CNP and ORCHIDEE-CNP. GOLUM-CNP stands for Global Observation-based Land-ecosystems Utilization Model of Carbon, Nitrogen and Phosphorus (GOLUM-CNP) v1.0, which is a data-driven model of steady-state C, N and P cycles for present-day (2001–2016) conditions. C, N and P losses via fire in ORCHIDEE-CNP are ignored. Numbers in square brackets indicate the standard deviations to account for the spatial spread of C, N and P fluxes.

4 Results

4.1 Carbon, nitrogen, and phosphorus flows and storages

We compared the simulated fluxes of C, N and P within natural ecosystems for the period 2001–2010 to the data-driven estimates from GOLUM-CNP (Table 1; Sect. S1B in Supplement) on the global scale and for natural ecosystems at biome scale. Modeled global C, N and P fluxes in ORCHIDEE-CNP are comparable with the estimates by GOLUM-CNP (Fig. 2). One exception is that ORCHIDEE-CNP simulates a fourfold lower P leaching from soils (3.7 ± 9.7 mg P m⁻² yr⁻¹) than GOLUM-CNP (23 mg P m⁻² yr⁻¹) (Fig. 2), which mainly occurs in forest ecosystems (Fig. S1). Note that GOLUM-CNP presents the steady-state C, N and P cycles in natural biomes, omitting human perturbations that have strongly altered the flows of C, N and P during the recent past. The impact of such perturbations on the nutrient flows are analyzed in detail in Sects. S2 and S3 in the Supplement.

In terms of C and nutrient storages, ORCHIDEE-CNP simulated comparable soil C, N and P storage (soil organic matter and litter) but higher vegetation C, N and P than GOLUM-CNP. Detailed comparisons for the spatial pattern of soil organic carbon (SOC) and forest aboveground C against observation-based datasets can be found in Figs. S2 and S3.

4.2 Resource use efficiencies

Here we evaluate the resource use efficiencies of GPP for light (L), water (W), C, N and P defined by

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\text{LUE}={\displaystyle \frac{\text{GPP}}{\text{f APAR}\times \text{PAR}}},\text{(3)}& {\displaystyle }& {\displaystyle }\text{WUE}={\displaystyle \frac{\text{GPP}}{\text{ET}}},\text{(4)}& {\displaystyle }& {\displaystyle }\text{CUE}={\displaystyle \frac{\text{GPP}}{\text{GPP}}},\text{(5)}& {\displaystyle }& {\displaystyle }\text{NUE}={\displaystyle \frac{\text{GPP}}{F_{\text{N}}}},\text{(6)}& {\displaystyle }& {\displaystyle }\text{PUE}={\displaystyle \frac{\text{GPP}}{F_{\text{P}}}},\end{array}$$

where GPP is the annual gross primary productivity (g C m⁻² yr⁻¹), fAPAR the fraction of absorbed photosynthetically active radiation (%), PAR the annual photosynthetically active radiation (W m⁻² yr⁻¹), ET the annual evapotranspiration (mm m⁻² yr⁻¹), and F_N and F_P the total N uptake (g N m⁻² yr⁻¹) and P uptake by plants (g P m⁻² yr⁻¹), respectively. We calculated fAPAR in ORCHIDEE-CNP and ORCHIDEE as a function of leaf area index (LAI): fAPAR $=\mathrm{1}-\exp (-\mathrm{0.5}\cdot \text{LAI})$ (Ito et al., 2004).

Compared to observed LUE (Sect. S1E in the Supplement), ORCHIDEE-CNP modeled median values at the biome level are generally lower but still within the ranges of uncertainties of observation-based datasets (Sect. 3.2) except for tropical (TRF) and temperate deciduous forest (TEDF). In comparison to ORCHIDEE, ORCHIDEE-CNP simulated LUEs closer to observations for four out of six biomes (TECF, BOCF, TEG, TRG) (Fig. 3a).

Compared to observed WUE, the ORCHIDEE and ORCHIDEE-CNP simulated values fall within the uncertainty range of observations (Fig. 3b). However, the WUE values from ORCHIDEE-CNP are on the high end of the range for temperate conifers (TECF) and BOCF and on the low end for temperate and tropical grasslands (TEG and TRG). The highest median WUEs were correctly simulated in temperate forests by ORCHIDEE-CNP (Fig. 3b), but the lowest WUE values were simulated in temperate instead of tropical forests.

Compared with observed CUE, ORCHIDEE-CNP simulated comparable values for TEDF and TECF but lower values for TRF, BOCF and grasslands. Both ORCHIDEE-CNP and ORCHIDEE cannot capture the increase in CUE from tropical to boreal forests apparent in the observation-based products (Fig. 3c) and in measurements from forest sites (Piao et al., 2010). In comparison to ORCHIDEE, ORCHIDEE-CNP simulated CUEs closer to observations for four out of six biomes (TEDF, TECF, BOCF, TEG) with respect to the median and spread.

Figure 3 Comparison of annual use efficiencies of light (LUE), water (WUE) and carbon (CUE) between ORCHIDEE-CNP, ORCHIDEE and satellite-based estimations for six biomes: tropical rainforest (TRF), temperate deciduous forest (TEDF), temperate conifer forest (TECF), boreal conifer forest (BOCF), temperate grass (TEG) and tropical grass (TRG). The whiskers indicate the interquartile (box) and 95 % confidence intervals (dashed lines). Grey boxes indicate the satellite-based estimations (referenced). The grey shaded areas indicate the uncertainties of resource use efficiencies given by referenced estimations, which involve uncertainties for multi-estimations and spatial variability for each estimation.

Consistent with site observations of NUE from Gill and Finzi (2016) and GOLUM-CNP outputs, ORCHIDEE-CNP correctly simulated the high values of TECF and the low values of tropical forests (Fig. 4a). However, compared with site observations of PUE from Gill and Finzi (2016), showing a PUE decrease from tropical to boreal regions, ORCHIDEE-CNP simulated a rather flat value (Fig. 4b). This suggests a P limitation that is too strong in high-latitude ecosystems, consistent with the fact that the model underestimates peak northern GPP and the northern land sink (Sect. 4.6). Nevertheless, the model-simulated PUE values fall in the range of GOLUM-CNP estimates. Tropical C₄ grasslands have higher simulated NUE and PUE than temperate C₃ grasslands, consistent with GOLUM-CNP (Fig. 4).

Figure 4 Violin plots of nitrogen use efficiency (NUE; a) and phosphorus use efficiency (PUE; b) by ORCHIDEE-CNP, GOLUM-CNP and observations (Gill and Finzi, 2018) for six biomes: tropical rainforest (TRF), temperate deciduous forest (TEDF), temperate conifer forest (TECF), boreal conifer forest (BOCF), temperate grass (TEG) and tropical grass (TRG). Open circles are the medians of all grid cells within each biome, with balloons representing the probability density distribution of each value. Black whiskers indicate the interquartile.

4.3 CO₂ fertilization effect

We compare the simulated response of plant productivity to increasing CO₂ during the historical period (i.e., CO₂ fertilization effect Eco₂) with observation-based estimates for C₃ plants from the historical change in deuterium isotopomers in leaf herbarium samples (Ehlers et al., 2015). For global (C₃ and C₄) vegetation we compare to indirect evidence from carbonyl sulfide (COS) atmospheric ice-core observations (Campbell et al., 2017). The CO₂ fertilization effect is defined here by the GPP ratio ($E_{{\mathrm{CO}}_{\mathrm{2}}}$):

$$\begin{array}{}\text{(7)}& E_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{{\text{GPP}}_{\mathrm{396}}}{{\text{GPP}}_{\mathrm{296}}}},\end{array}$$

where GPP₂₉₆ indicates pre-industrial GPP (g C m⁻² yr⁻¹) under a CO₂ concentration of 296 ppm and GPP₃₉₆ under a current CO₂ concentration of 396 ppm. Those CO₂ concentrations of 296 and 396 ppm correspond to tropospheric mixing ratios of CO₂ in the years ∼ 1900 and 2013, respectively, similar to values used for estimating the response of GPP to a ∼ 100 ppm CO₂ increase in Ehlers et al. (2015) and Campbell et al. (2017).

Figure 5 Comparisons between pre-industrial GPP with an atmospheric CO₂ concentration of 296 ppm (GPP₂₉₆) and current GPP with an atmospheric CO₂ concentration of 396 ppm (GPP₃₉₆) for all natural plants (a, b) and natural C₃ plants (c, d) by ORCHIDEE-CNP (a, c) and ORCHIDEE (b, d). The color scale shows the point density. Different point densities and patch sizes for ORCHIDEE and ORCHIDEE-CNP are due to the different spatial resolution (2^∘ × 2^∘ for ORCHIDEE-CNP and 0.5^∘ × 0.5^∘ for ORCHIDEE). The ratio between GPP₃₉₆ and GPP₂₉₆ indicates the CO₂ fertilization effects ($E_{{\mathrm{CO}}_{\mathrm{2}}}$). Green dashed areas indicate the observed $E_{{\mathrm{CO}}_{\mathrm{2}}}$ from the Campbell et al. (2017) COS records. Pink lines indicate the observed $E_{{\mathrm{CO}}_{\mathrm{2}}}$ from Ehlers et al. (2015).

Figure 6 Global pattern of N (O_N, a, b) and P openness and that (O_P, c, d) simulated by ORCHIDEE-CNP (a, c) and GOLUM-CNP (b, d). Pixels with managed lands > 50 % in ORCHIDEE-CNP were masked. The same area was masked from the pattern of O_N and O_P for GOLUM-CNP.

Modeled $E_{{\mathrm{CO}}_{\mathrm{2}}}$ by ORCHIDEE-CNP of natural biomes ranges between 1.0 and 1.3 for most regions (Fig. 5a), slightly lower than global $E_{{\mathrm{CO}}_{\mathrm{2}}}$ derived from COS of 1.26–1.36 (Campbell et al., 2017). Modeled $E_{{\mathrm{CO}}_{\mathrm{2}}}$ values for C₃ plants (Figs. 5c, S4) are also consistent with $E_{{\mathrm{CO}}_{\mathrm{2}}}$ from herbarium samples (Ehlers et al., 2015), equal to 1.23. When compared to ORCHIDEE without nutrient cycles, we found that ORCHIDEE-CNP simulates smaller and more realistic values of $E_{{\mathrm{CO}}_{\mathrm{2}}}$ (Fig. 5c, d) but with lower values in boreal regions that could not be checked against observations (Fig. S5).

4.4 Ecosystem nutrient openness and nutrient turnover

Nutrients taken up by plants are either recycled within the ecosystem or acquired from external sources (P weathering of primary and secondary minerals, atmospheric N and P deposition, biological nitrogen fixation – BNF, and N and P fertilizer addition to cultivated lands). Wang et al. (2018) calculated an indicator of the openness of N and P cycling in natural ecosystems as the ratio of external inputs of N and P into the ecosystem to the total amount of N and P that plants use for GPP. Similarly, we diagnosed the openness for N and P (O_N and O_P) from the ORCHIDEE-CNP output by

$$\begin{array}{}\text{(8)}& O_x={\displaystyle \frac{I_x}{F_x+{\text{RSB}}_x}},\end{array}$$

where I_x is the annual external nutrient input (g X m⁻² yr⁻¹), F_x the annual plant uptake of soil nutrients (g X m⁻² yr⁻¹) and RSB_x the flux of nutrients recycled within plants (g X m⁻² yr⁻¹) by foliar nutrient resorption prior to leaf shedding. External nutrient inputs include atmospheric N deposition and BNF, as well as P deposition and P release from rock weathering.

Modeled O_N in natural biomes by ORCHIDEE-CNP showed only a small variance across the globe, whereas GOLUM-CNP predicts a higher O_N in tropical and temperate regions than in boreal regions (Fig. 6a, b). O_P values are below 15 % in ORCHIDEE-CNP for most biomes, of a similar order of magnitude as in GOLUM-CNP (Fig. 6c, d). ORCHIDEE-CNP simulates a lower O_N in tropical natural biomes than GOLUM-CNP, which is mainly due to lower but more realistic tropical BNF in ORCHIDEE-CNP compared to GOLUM-CNP (Sect. S4 in the Supplement). ORCHIDEE-CNP simulates a higher O_N in high-latitude grassland (Fig. 6a, b) than GOLUM-CNP, which is due to overestimation of BNF in NH in ORCHIDEE-CNP (Sect. S4 in the Supplement). Modeled O_P in natural biomes by ORCHIDEE-CNP compares well with GOLUM-CNP except for central Africa (Fig. 6c, d). This is primarily because ORCHIDEE-CNP used a lower P deposition forcing than GOLUM-CNP.

Residence time quantifies the average time it takes for an N (or P) molecule to enter and leave the ecosystem (τ_N and τ_P). In this study, we adopted the approach of Carvalhais et al. (2014) for the carbon residence time. We define the residence time of N and P as the ratio of total respective nutrient stock in the ecosystem to their respective total input flux:

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\text{N}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^{\mathrm{5}}{\text{N}}_i+{\text{N}}_{\text{inorg}}}{{\text{N}}_{\text{d}}+\text{BNF}}},\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\text{P}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^{\mathrm{5}}{\text{P}}_i+{\text{P}}_{\text{inorg}}}{{\text{P}}_{\text{d}}+{\text{P}}_{\text{w}}}},\end{array}$$

where N_i indicates the N mass (g N m⁻²) in organic matter pools i (with i= plant, litter, SOM pools), N_inorg is the sum of all inorganic N pools, and N_d and BNF are N deposition and biological N fixation rates, respectively (g N m⁻² yr⁻¹). Similarly, P_i is the P mass (g P m⁻²) in organic matter pools, P_inorg the sum of inorganic P pools, and P_d and P_w are P deposition and P weathering release rates, respectively (g P m⁻² yr⁻¹).

Figure 7 Violin plots of the residence time of (a) N and (b) P cycles for six biomes: tropical rainforest (TRF), temperate deciduous forest (TEDF), temperate coniferous forest (TECF), boreal coniferous forest (BOCF), temperate grass (TEG) and tropical grass (TRG). Open circles are the medians of all grid cells within each biome, with balloons representing the probability density distribution of each value. Black whiskers indicate the interquartile.

Modeled median τ_N of natural biomes in ORCHIDEE-CNP varies between 56 and 1585 years, while τ_P varies within a large range of 101 to 223 870 years (Fig. 7). ORCHIDEE-CNP captured the order of magnitude of τ_N and τ_P for forests found in GOLUM-CNP. Longer median τ_N (1585 years) and τ_P (1 223 870 years) are simulated for boreal forest compared to temperate and tropical forests (251–794 years for τ_N and 891–7080 years for τ_P) and grassland (56–158 years for τ_N and 101–468 years for τ_P) by ORCHIDEE-CNP, consistent with results from GOLUM-CNP. However, for grasslands, simulated τ_N (56-158 years) and τ_P (101–468 years) are 5–11-fold shorter than in GOLUM-CNP (Fig. 7).

Figure 8 Comparisons of the leaf N : P ratio between ORCHIDEE-CNP, data-driven estimates and observations. (a) The global pattern of the mean leaf N : P ratio over 2001–2016 for ORCHIDEE and (b) for mean leaf N : P in Butler et al. (2017). Panels (c) and (d) are the 25th and 75th percentile, respectively, of the leaf N : P ratio by Butler et al. (2017). Dots in (a) indicate the area with a leaf N : P ratio from ORCHIDEE-CNP falling into the 25th–75th percentiles of the Butler et al. (2017) estimation. (d) The latitude distributions of the leaf N : P ratio for ORCHIDEE-CNP, the Butler et al. (2017) estimation and site measurements. Red shared area indicates the uncertainty from latitudinal spreads of the leaf N : P ratio for ORCHIDEE-CNP. The grey shaded area indicates the uncertainty from both the estimations and latitudinal spreads for Butler et al. (2017). Blue and yellow lines indicate the 10th and 90th percentiles, respectively, of measured leaf N : P ratios in each bin of 3^∘ latitude.

4.5 Stoichiometry

4.5.1 Foliar stoichiometry

Leaf N : P ratios for natural biomes predicted by ORCHIDEE-CNP vary between 15 and 25 (Fig. 8a). The observed decline in median leaf N : P ratios with increasing latitude was not reproduced by the model (Sect. S1E1 in the Supplement; Fig. 8e), although the modeled latitudinal distribution of leaf N : P ratios remained within the 10–90th quantiles of the site-level data (Kerkhoff et al., 2005; McGroddy et al., 2004; Reich and Oleksyn, 2004). Further, the simulated leaf N : P ratios fall within the interquartile of upscaled site measurements by Butler et al. (2017) for most of the globe, with the exception of regions north of 55^∘ N where leaf N : P ratios are outside the observation-based range, suggesting a P constraint relative to N that is too strong (Fig. 8).

4.5.2 Soil stoichiometry

Here we evaluate the modeled C : N, C : P and N : P ratios of soil organic matter for different biomes against data from the large compilation of measurements for soils (0–60 cm depth) by Tipping et al. (2016). Modeled C : N ratios fall into much more narrow ranges (7.8–11.8 for the widest interquartile range) compared to the observations (11.1–20.5; Fig. 9a) as a result of prescribing constant C : N ratios in ORCHIDEE-CNP (Goll et al., 2017a). SOM P content varies in ORCHIDEE-CNP as a consequence of varying biochemical phosphorus mineralization rates (Sect. S7), and thus C : P and N : P ratios of SOM show pronounced variation in space. ORCHIDEE-CNP simulates comparable N : P ratios as measurements in terms of both the median value and distributions for tropical forests, but it overestimates the observed N : P ratios by 108 %–327 % in temperate forest, tropical and temperate grassland soils (Fig. 9b, c). The higher observed C : P and N : P in forest compared to grassland soils are not captured by ORCHIDEE-CNP (Fig. 9b, c). We also compared ORCHIDEE-CNP N : P ratios to the results of GOLUM-CNP, which were based on data from Zechmeister-Boltenstern et al. (2015) that are more limited than Tipping et al. (2016), and found an overestimation for temperate forests, tropical forests and temperate grasslands.

Figure 9 C : N, C : P and N : P ratios of soil organic matter by ORCHIDEE-CNP and plot-level measurements by Tipping et al. (2016) for four biomes: tropical forest (TRF), temperate forest (TEF), tropical grass (TRG) and temperate grass (TEG). Soil C : N : P ratios for ORCHIDEE-CNP are calculated for the total soil pool and include soil passive, slow and active pools, while measurements by Tipping et al. (2016) are for soils of 0–60 cm depth. The letters “a”, “b” and “c” indicate the significance of differences among biomes from the analysis of variance (ANOVA).

4.6 Nutrient effects on carbon cycling

We analyze the performance of ORCHIDEE-CNP v1.2 and ORCHIDEE without nutrient cycles with respect to the spatiotemporal patterns of GPP, NPP and net biome productivity.

Global GPP and NPP simulated by ORCHIDEE-CNP averaged over the period 2001–2010 are 119 and 48 Pg C yr⁻¹, respectively, which are both within ranges of the data-driven products listed in Table 1 (Sect. S1C in the Supplement; Table S2). GPP and NPP simulated by ORCHIDEE-CNP are lower than those simulated by ORCHIDEE (140 Pg C yr⁻¹ for GPP and 60 Pg C yr⁻¹ for NPP). The values from ORCHIDEE are on the high end of the range of estimates from the data-driven products in Table 1. ORCHIDEE-CNP simulated comparable GPP values for most parts of the globe (Fig. S6a) and comparable NPP values for most of the northern high latitudes (Fig. S6b), which lie within the range given by the data-driven products.

Figure 10 The performances of ORCHIDEE and ORCHIDEE-CNP on the interannual variability of de-trended anomalies of GPP during 2001–2010 (a), the seasonal variability of mean GPP across 2001–2010 (b), and the interannual variability of net biome productivity (NBP) (c). Two statistics were used to represent the model performance: coefficient of determination (R²) and relative mean square deviation (rMSE). For (a) and (b), the evaluations are for the globe, the Northern Hemisphere (30–90^∘ N; NH), the northern tropics (0–30^∘ N; NT), the southern tropics (0–30^∘ S; ST) and the Southern Hemisphere (30–90^∘ S; SH). Two sets of observation-based GPP products (BESS GPP and MTE GPP) were used for the comparison. For (c), the evaluations are for the globe, the Northern Hemisphere (30–90^∘ N, NH), the tropics (30^∘ S–30^∘ N) and the Southern Hemisphere (30–90^∘ S; SH). The mean value across TRENDY ensemble models (v6) and three sets of NBP from inversion datasets were used as the reference database for the comparison with different available periods (TRENDY Ensemble: 1959–2016; CAMS: 1979–2016; JENA: 1985–2016; CTracker: 2001–2016).

Interannual and seasonal variations of GPP reflect the response of ecosystems to interannual or seasonal climatic variability, as well as the effects of natural (e.g., fires, wind throw, insect outbreaks and storms) and anthropogenic disturbances (e.g., land management and land cover change) (Anav et al., 2015). Regarding the interannual anomalies of de-trended GPP (GPP_int) for the period 2001–2011, estimations on a global scale from ORCHIDEE-CNP show rather good correlation with the observation-driven model BESS GPP (R²=0.71) but not with MTE GPP (R²=0.11) (Fig. 10a). ORCHIDEE performs somewhat worse on a global scale than ORCHIDEE-CNP, primarily due to its low performance in the NH. We find that inclusion of nutrients in ORCHIDEE leads to a lower model prediction error on a global scale and for all latitudinal bands irrespective of the observation-based product (Fig. 10a).

Regarding the seasonal variation of GPP over the period 2001–2011, the predictions of ORCHIDEE-CNP are in good agreement with observation-based estimates and show no significant differences when compared to ORCHIDEE, except for tropical regions (Fig. 10b). Here, the model errors in seasonal variations of GPP are substantially larger for ORCHIDEE-CNP than for ORCHIDEE (Fig. 10b).

Net biome productivity (NBP) is defined as the net C exchange between the atmosphere and the terrestrial biosphere, which is the sum of net primary productivity, heterotrophic respiration and emissions due to disturbances; positive values denote a land carbon sink. Compared to the three sets of atmospheric inversions (CAMS, JENA and CTracker), ORCHIDEE(-CNP) performs slightly worse than the mean of predictions from 16 land surface models from TRENDY ensembles (v6) (Fig. 10c). ORCHIDEE-CNP shows a worse performance in interannual variability of NBP than ORCHIDEE when compared against inversion datasets at global scale and for the Northern Hemisphere. However, ORCHIDEE-CNP improved the performance of the interannual variability of NBP against inversion datasets relative to ORCHIDEE for tropical regions (higher R² and lower rMSE), with closer or even better fitness against inversion datasets than the mean value of TRENDY ensemble models (Fig. 10c).

5 Discussion

We performed a detailed evaluation of ORCHIDEE-CNP in terms of four nutrient-related ecosystem properties that control ecosystem gas exchanges and carbon storage: vegetation resource use efficiencies, CO₂ fertilization effect, ecosystem N and P turnover and openness, and large-scale pattern of ecosystem stoichiometries.

We find that the inclusion of nutrients tends to lead to improvements in simulated resource use efficiency of plant resources (light, carbon, water) on a biome scale (Sect. 5.1). In line with changes in resource use efficiency, the sensitivity of GPP to variations in climate is improved, leading to improved interannual variation in GPP, in particular for the Northern Hemisphere (Sect. 5.5). In addition, the response of GPP to an increasing atmospheric CO₂ concentration is improved (Sect. 5.2). However, model biases in C fluxes remained or increased, for example, in the NBP of the Northern Hemisphere. The analysis of nutrient use efficiencies (Sect. 5.1), stoichiometry (Sect. 5.4), and ecosystem openness and turnover of nutrients (Sect. 5.3) reveals biases in boreal regions that might be related to issues with soil organic matter accumulation that is too strong and the dependency of photosynthesis on leaf nutrients in needleleaf PFTs. On a seasonal scale, we found a general deterioration of the simulated seasonal cycle of GPP due to the inclusion of nutrient cycles (Sect. 5.5).

In the following, we discuss in more detail the model performance with respect to nutrient cycles and their effects on simulated C fluxes, and we propose ways to address model biases.

5.1 Inclusion of nutrient cycling improves use efficiencies of other plant resources

Resource use efficiency (RUE) is an ecological concept that measures the proportion of supplied resources that support plant productivity; i.e., it relates realized to potential productivity (Hadapp et al., 2019). It is therefore a critical ecosystem property that relates resource availability to ecosystem productivity, as well as being affected by resource availability.

With the inclusion of the additional plant resources nitrogen and phosphorus, changes in the simulated vegetation use efficiencies of resources like water (WUE), light (LUE) and carbon (CUE) are expected. Indeed, the annual use efficiencies on the biome scale differ between ORCHIDEE-CNP and ORCHIDEE. In comparison to observation-based estimates, the inclusion of nutrient cycles tends to improve simulated LUE and CUE and WUE (Fig. 3).

Both ORCHIDEE-CNP and ORCHIDEE generally underestimate annual LUE for forest biomes (Fig. 3a), which is due to a high bias in fAPAR in both models (28 %–380 % for ORCHIDEE and 80 %–173 % for ORCHIDEE-CNP) (Fig. S4a, b). Although the bias in LUE for TRF is higher, the bias in GPP is largely reduced, whereas the bias in fAPAR is similar in ORCHIDEE-CNP compared to ORCHIDEE (Fig. S4a, b), indicating general issues in ORCHIDEE with respect to how light is transferred within the canopy in tropical forest. Both versions assume a constant canopy light extinction coefficient of 0.5, omitting variations among biomes due their distinctive canopy architectures (Ito et al., 2004). Improving this part of the model requires a canopy light transfer scheme that better accounts for canopy structure (Naudts et al., 2015) and the inclusion of different light components including diffuse incoming, scattered and direct light (Zhang et al., 2020).

ORCHIDEE-CNP simulated a lower WUE than ORCHIDEE with the exception of coniferous biomes (Fig. 3b). The improvement of WUE in TRF is related to improvements in GPP and ET, while the overestimation of WUE in coniferous-dominated biomes by ORCHIDEE-CNP is related to an overestimation of GPP (Fig. S4c). The latter is likely a result of the application of a relationship between photosynthetic capacity and the leaf nutrient concentration, which is based on measurements from broadleaf species for all PFTs. Kattge et al. (2009) showed that coniferous PFTs have a ∼ 40 % lower carboxylation capacity for a given leaf nitrogen concentration than other PFTs. The omission of this could explain the bias in coniferous GPP in ORCHIDEE-CNP. Uncertainties in evaluation datasets hamper a more detailed evaluation of the variations of WUE among biome types.

We found that the inclusion of nutrient cycles improved the spatial variability in simulated CUE, but general biases remain (Fig. 3c), and uncertainties in observation-based estimates are large. Improvements are mainly found in temperate biomes (TEDF, TECF and TEG), indicating that the allocation of GPP to respiration and biomass growth, which is controlled by nutrient availability, works reasonably well. ORCHIDEE-CNP underestimates CUE for tropical biomes (TRF and TRG) more strongly than ORCHIDEE, despite substantially reduced biases in NPP and GPP (Fig. S4d). However, we should be cautious in drawing conclusions considering the large uncertainty in MODIS CUE (He et al., 2018).

NUE and PUE on the biome scale compare well to estimates (Fig. 4), indicating that ORCHIDEE-CNP is able to simulate the coupling strength between C, N and P cycles. However, ORCHIDEE-CNP underestimates PUE in tropical forests. A sensitivity analysis by GOLUM-CNP indicated that NUE and PUE were most sensitive to the NPP allocation fractions (especially to woody biomass) and foliar stoichiometry (Wang et al., 2018). Therefore, we attribute the biases in PUE to the biases in foliar stoichiometry (Fig. 8) and to issues in plant internal P allocation in ORCHIDEE-CNP (Fig. S1).

5.2 Inclusion of nutrient cycling improves CO₂ fertilization effect

The effect of CO₂ fertilization on terrestrial ecosystem productivity is thought to be the dominant driver behind the current land carbon sink. The strength of the fertilization effect on GPP differs strongly between LSMs (Friedlingstein et al., 2014). We used proxies of the historical increase in GPP for an indirect model evaluation of the CO₂ fertilization effect from COS and deuterium measurements of herbarium samples (Ehlers et al., 2015; Campbell et al., 2017), and we found that ORCHIDEE-CNP has smaller and more realistic $E_{{\mathrm{CO}}_{\mathrm{2}}}$ than the same model without nutrients (Fig. 5), in particular for C₃ plants and in boreal regions (Fig. S5). Both ORCHIDEE-CNP and ORCHIDEE simulated an $E_{{\mathrm{CO}}_{\mathrm{2}}}$ for C₄ grass of ∼ 1, as the carboxylation of C₄ plants is weakly influenced by elevated CO₂ (Osmond et al., 1982; Pearcy and Ehleringer, 1984; Bowes, 1993). This indicates that the inclusion of N and P constraints on GPP leads to a more realistic CO₂ fertilization effect in ORCHIDEE-CNP.

5.3 Ecosystem nutrient turnover and openness indicates model biases in boreal phosphorus availability

The capacity of ecosystems to sequester and store additional carbon depends on their ability to supply nutrients for the buildup of organic matter. Enhanced internal nutrient recycling and the accumulation of nutrients over time in ecosystems are theoretically possible mechanisms through which nutrients can be supplied. Therefore, it is important for simulating changes in land carbon storage on decadal timescales and longer that models capture the dependency of ecosystem production to external nutrient sources (i.e., openness of N and P cycles) (Cleveland et al., 2013) and the residence time of nutrients within ecosystems. Besides being related to each other, openness and residence times are also related to the inflows and outflows of nutrients (Eqs. 9 and 10) as well as the turnover time of nutrients in specific ecosystem compartments.

We find that ORCHIDEE-CNP simulates the openness of nutrient cycles, including differences among biomes that are close to estimates from the model–data fusion framework GOLUM-CNP (Fig. 6; Sect. 4.4). There are differences in the openness of N (O_N) in tropical natural biomes and the openness of P (O_P) in central Africa, which are related to lower but more realistic tropical BNF in ORCHIDEE-CNP (Sect. S4) and a difference in the prescribed P deposition compared to GOLUM-CNP. Simulated nutrient losses due to aquatic transport are generally in good agreement with independent estimates (Sect. S5).

Residence times of N and P (τ_N and τ_P) in ORCHIDEE-CNP compare generally well to estimates from GOLUM-CNP: ORCHIDEE-CNP simulates shorter τ_N and τ_P in tropical and temperate biomes compared to boreal ones, in line with GOLUM-CNP (Fig. 7). This indicates that ORCHIDEE-CNP is able to reproduce large-scale patterns in the nutrient residence time of biomes, with one exception. In boreal regions, we find that ORCHIDEE-CNP simulates higher τ_P for BOCF due to the higher standing P stocks of biomass and soil organic matter than in GOLUM-CNP (Fig. S1). This indicates that ORCHIDEE-CNP is likely underestimating P availability in boreal regions. The underlying processes of biochemical P mineralization (Sect. S7) and sorption of P to soil particles (Sect. S6) are reasonably well captured in ORCHIDEE-CNP.

5.4 Model biases in stoichiometry indicate need for refinement of process representation

Leaf and soil stoichiometry are key indexes to characterize the relative ecosystem N and P limitation (e.g., Güsewell, 2004). Measurements show a decrease in foliar N : P ratios from low to high latitudes in natural ecosystems (McGroddy et al., 2004; Reich and Oleksyn, 2004; Kerkhoff et al., 2005). This is seen as evidence for tropical vegetation being generally more P- than N-limited, in contrast to extratropical vegetation (Reich and Oleksyn, 2004). The observed trend of foliar N : P ratios was not reproduced by ORCHIDEE-CNP (Fig. 8), which simulated a flat foliar N : P latitudinal profile. In contrast to the majority of global models, wherein leaf N : P ratios are either prescribed (Goll et al., 2012) or vary within a PFT-specific range (Wang et al., 2010), we conservatively assumed a globally uniform range to let the model freely calculate leaf N : P stoichiometry. It is not trivial to pin down the failure of the model to capture the latitudinal trend in leaf N : P ratios, which could be due to (1) omitted variability in leaf P resorption efficiencies, which varies among biomes between 46 % and 66.6 % (Reed et al., 2012) but was set to 65 % in ORCHIDEE-CNP, (2) the simplistic parameterization of nutrient investment into different plant tissues, (3) and the omission of the diversity of nutrient acquisition pathways (e.g., mycorrhizal association) and rooting strategies (Warren et al., 2015). Testing new formulations for plant growth based on optimality principles (Kvakić et al., 2020) and the refinement of nutrient acquisition pathways (Sulman et al., 2017) are ways forward to improve the model.

Regarding soil stoichiometry, measurements show that tropical biomes have lower soil C : N and higher soil C : P and soil N : P than temperate biomes (Tipping et al., 2016), echoing the pattern of leaf stoichiometry. ORCHIDEE-CNP fails to capture these patterns (Fig. 9). Modeled soil N : P and C : P for tropical forests are comparable to measurements but are too low in temperate forest, tropical forest and temperate grass, which is most likely related to a nutrient immobilization that is too strong in accumulating soil organic matter (Figs. S1); this tends to push systems into P limitation rather than N limitation as O_N is larger than O_P (Fig. 6). In general, the spread in soil P concentration is well represented by ORCHIDEE-CNP. The rudimentary representation of organic matter decomposition and the lack of nutrient effects on decomposer carbon use efficiency (see Zhang et al., 2018, for possible improvements; Sect. 5.5) are likely contributing to the biases. New developments, including explicit representation of decomposer communities and soil organic matter stabilization (Zhang et al., 2020), will be included in the next model version.

5.5 Nutrient effects on carbon cycling

In the following we discuss the implications for the simulated carbon fluxes of changes in plant resource use efficiencies and the sensitivity of plant productivity to increasing CO₂ due to the inclusion of nutrient cycles. We link biases in the simulated carbon fluxes to biases in nutrient cycling, which allows us to prioritize follow-up model development.

5.5.1 Inclusion of nutrient cycling improves the interannual variability of GPP

To what extent nutrient effects on vegetation affect the sensitivity of ecosystem CO₂ fluxes to climatic variation is unclear (Goll et al., 2018). For instance, drought can reduce nutrient use by decreasing GPP, but it also slows down decomposition, which supplies nutrients for plant uptake. Further, N : P stoichiometry is also strongly modified by drought and warming towards increased N : P in whole plant biomass (Yuan and Chen, 2015). Here we found that the inclusion of N and P cycles in ORCHIDEE affects the interannual variability of GPP for all vegetation types. In ORCHIDEE-CNP, the interannual variation (IAV) of GPP is better correlated with that of observation-based datasets than in ORCHIDEE globally and for the NH, but less correlated for other regions (Fig. 10a). Observation-based GPP estimates are uncertain, as some of them ignore soil-moisture-induced reductions of GPP during drought (Stocker et al., 2019), as well as soil thaw and snow-related effects (Jiang and Ryu, 2016). Thus, at the moment, it is difficult to falsify one model version over another and to constrain nutrient effects on the variation of GPP based on current observation-based GPP.

In order to further explore the underlying reasons for the general improvement in the IAV of GPP due to the inclusion of nutrient cycles, we analyzed the sensitivity of GPP anomalies to anomalies of temperature (S_T), precipitation (S_P) and shortwave radiation (S_R), all with mean annual values (Sect. S1I). We found that S_P by ORCHIDEE-CNP compares well with BESS GPP and MTE GPP, while it is overestimated in ORCHIDEE (Figs. S7 and S8). Thus, the difference in S_P is likely the major reason for the differences in IAV in NH between model versions, as S_T and S_R show only minor differences there. This provides confidence that the improvement in IAV of GPP in the NH is due to an improved sensitivity towards a climatic driver (i.e., S_P). For tropical regions, ORCHIDEE-CNP simulates more realistic S_P but higher biases in S_R than in ORCHIDEE, while observation-based estimates of S_T disagree on the sign and model versions show only minor differences (Fig. S7). Therefore, the deterioration of the IAV of tropical GPP by the inclusion of nutrient cycles is likely caused by enhanced biases in S_R due to a lowering of LUE of GPP (Sects. 4.2 and 5.1).

5.5.2 Inclusion of nutrient cycling deteriorates phenology and seasonality of GPP

The performance in reproducing seasonal variations of GPP was deteriorated by the inclusion of N and P nutrient cycles in ORCHIDEE (Fig. 10b). We found that biases in GPP are related to biases in the seasonality of the LAI introduced in ORCHIDEE-CNP (Figs. S9a and S10a). For the NH, the delayed increase in LAI in ORCHIDEE-CNP could be partly caused by nutrient shortage during the first half of the growing season, as indicated by the increasing leaf nutrient concentration throughout the growing season (Fig. S11). Several factors could lead to a supply of nutrients that is too low at the beginning of the growing season: an insufficient internal plant nutrient reserve due to resorption of nutrients that is too low prior to leaf shedding or an underestimation of nutrient uptake during the dormant season, an insufficient investment into root growth to acquire nutrients, and an overestimation of soil nutrient losses during the dormant season leaving the soil nutrient depleted at the beginning of the growing season. Many of the related processes (e.g., root phenology, mineralization, nutrient resorption, growth allocation) are only rudimentarily represented. For tropical regions, ORCHIDEE-CNP simulates a quasi-flat seasonal cycle of GPP, in contrast to a peak of GPP during the wet season in MTE GPP and BESS GPP, which is correctly captured by ORCHIDEE (Fig. S9b, c). The reduction of seasonal GPP in ORCHIDEE-CNP compared to ORCHIDEE is more pronounced in the dry season (∼ 100 g C m⁻²) than in the wet season (Fig. S9b, c), concurrent with a larger reduction of LAI in the dry season (Fig. S10b, c). Tropical phenology is currently only rudimentarily represented in ORCHIDEE(-CNP) (Chen et al., 2020), causing a suboptimal allocation of nutrients to leaves that could cause the biases in the seasonal cycle of GPP and LAI. Model–data assimilation of phenology (Williams et al., 2009; MacBean et al., 2018; Bacour et al., 2019) and efforts to better characterize processes related to plant resource investment into different tissues and symbioses (Prentice et al., 2015; Warren et al., 2015; Jiang et al., 2019) as well as leaf age effects during the year for evergreen forests (Chen et al., 2020) should help to reduce tropical phenology biases in future versions of ORCHIDEE-CNP.

5.5.3 Inclusion of nutrient cycling leads to an underestimation of the land carbon sink

Current LSMs unanimously conclude that CO₂ fertilization is the main driver of the land carbon sink and its trend (Friedlingstein et al., 2014), but it remains unclear to what extent other drivers (i.e., climate change, land management, nutrient deposition) contribute to the sink as well. Also, it remains unclear how commonly omitted dynamics (climate and management induced effects on tree mortality, nutrients) lead to overestimation of the contribution of CO₂ fertilization in models (Ellsworth et al., 2017; Fleischer et al., 2019). ORCHIDEE-CNP simulates a land carbon sink over the past decades that is lower than other dynamic global vegetation models (DGVMs) and atmospheric inversions (Fig. S12), despite the fact that the response of GPP to CO₂ in ORCHIDEE-CNP is in line with proxy data (Fig. 5; Sect. 5.2). In particular, the NH carbon sink, which has persistently increased over the last 50 years (Ciais et al., 2019), is strongly underestimated. The few Free Air Carbon Enrichment (FACE) studies that have experimentally applied elevated CO₂ levels in mature stands found no increase in biomass production (Bader et al., 2013; Klein et al., 2016; Körner et al., 2005; Sigurdsson et al., 2013; Ellsworth et al., 2017); thus, an increase in GPP does not necessarily translate into an increase in biomass production, whereas in most DGVMs wherein mortality is constant and growth follows GPP, biomass production is inevitably coupled to GPP. Based on upscaling of data from FACE experiments, Terrer et al. (2019) suggested that the effect of elevated CO₂ on biomass may be severely overestimated (on average by a factor of 3.6) in LSMs that ignore nutrients. It would be tempting to conclude from this study that ORCHIDEE-CNP is “right” in its underestimation of the carbon sink, whereas other models are “wrong” because they miss processes such as forest regrowth (Pugh et al., 2019) from, e.g., decreased harvesting pressure (Ciais et al., 2008) and thus have a realistic NH land sink for the wrong reasons. We also showed that ORCHIDEE-CNP underestimates peak GPP (Fig. S12b) and overestimates P limitations in the NH (Sect. 5.1, 5.3 and 5.4); thus, another explanation is that the NH sink in this study is too low because of P limitations that are too strong in this region. These two hypotheses explaining why we underestimate the NH sink (missing forest regrowth vs. overly strong nutrient limitations in the NH) are examined below.

The overly small NH carbon sink in ORCHIDEE-CNP may be explained by an immobilization of nutrients that is too strong in accumulating nutrient-rich organic matter, which leads to a reduction of plant-available nutrients, the so-called “progressive nutrient limitation” proposed by Luo et al. (2004), and subsequently to reduced biomass production. The amount of accumulated N and P immobilized into SOM in the NH during 1850–2016 reaches up to 75.3 g N m⁻² and 2.4 g P m⁻², respectively, which is twice as much as the accumulated respective nutrient inputs to ecosystems in this region during the same period (37.8 g N m⁻² and 1.6 g P m⁻²; Figs. S13 and S14). This suggests a strong progressive nutrient limitation in the model. The omission of nutrient controls on litter and SOM decomposition in the soil module of ORCHIDEE-CNP could have favored the immobilization of nutrients in accumulating SOM (Zhang et al., 2018). Microbe incubation and N fertilization experiments showed that a low availability of nutrients can hamper the buildup of SOM as more carbon gets respired by decomposers due to elevated energetic requirements of processing low-quality substrate (Recous et al., 1995; Janssens et al., 2010; Allison et al., 2009) and overall lower microbial activity (Wang et al., 2011; Knorr et al., 2005). Uncertainties with respect to the capability of ecosystems to up-regulate P mineralization when P becomes scarce could have contributed to the decline in plant-available nutrients with increasing SOM stocks. The inclusion of nutrient effects on decomposition and microbial dynamics in ORCHIDEE-CNP is ongoing (Zhang et al., 2018, 2020), but the lack of quantification of the ability of ecosystems to enhance P recycling hampers model developments.

The overly small NH carbon sink in ORCHIDEE-CNP may also be explained by the lack of representation of effects of forest age and management on biomass turnover and biomass production efficiency (i.e., CUE). Pugh et al. (2019) found that old-growth forests in the NH have a much smaller C sink than regrowing forests (< 0.1 Pg C yr⁻¹ compared to 0.86 Pg C yr⁻¹) for the period 2001–2010. Forest management effects on biomass production efficiency and biomass turnover are only rudimentarily represented in ORCHIDEE(-CNP). ORCHIDEE-CNP prescribes constant tree mortality rates (i.e., the fraction of total carbon in wood lost to litter), whereas in reality tree mortality rates change with management and climate conditions (Peng et al., 2011). Moreover, ORCHIDEE(-CNP) omits the effect of forest age on C uptake. Compared to data-driven estimates for C storage (Sect. S1G and S1H), ORCHIDEE-CNP simulates a higher global aboveground forest biomass (387 Pg C; 283 Pg C for GlobBiomass, Santoro, 2018; 221 Pg C for GEOCARBON, Operational Global Carbon Observing System, Avitabile et al., 2016; Fig. S2) but lower global soil organic carbon (801 Pg C; 4387 Pg C for Soilgrids, Hengl et al., 2017; 1680 Pg C for GSDE – Global Soil Dataset for use in Earth system models, Shangguan et al., 2014; Fig. S3).

6 Concluding remarks

In this study, we evaluated the performance of ORCHIDEE-CNP and found that the model has sufficient skills in capturing observed patterns in (1) vegetation resource use efficiencies, (2) CO₂ vegetation fertilization, (3) ecosystem N and P openness and turnover, and (4) leaf and soil stoichiometry. The inclusion of nutrients improves the simulation of the sensitivity of plant productivity to increasing CO₂ and to interannual variation in precipitation. However, the nutrient-enabled version cannot capture the current land carbon sink in the NH. This suggests that either the land carbon sink might be less a consequence of the CO₂ fertilization effect than of other processes that are currently not well resolved in global models (e.g., biomass turnover, land management) or that ORCHIDEE-CNP underestimates the ability of ecosystems to circumvent nutrient constraints on biomass built up under elevated CO₂. We propose the following focus to improve ORCHIDEE in the next model versions: (1) refine the canopy light absorption processes; (2) use model–data assimilation frameworks (like ORCHIDAS) to better calibrate root phenology, mineralization, nutrient resorption and growth allocation; (3) better represent soil processes related to decomposition, stabilization of soil organic matter (e.g., Zhang et al., 2018, 2020) and inorganic P transformation (e.g., Helfenstein et al., 2020); and (4) refine the dynamics of biomass turnover and biomass production efficiency, including effects of forest management and climate. Continued improvements of nutrient cycle representations will further reduce uncertainties in predicting the land carbon sink under climate change and rising atmospheric CO₂.