A new model of the coupled carbon, nitrogen, and phosphorus cycles in the terrestrial biosphere (QUINCY v1.0; revision 1996)

The dynamics of terrestrial ecosystems are shaped by the coupled cycles of carbon, nitrogen, and phosphorus, and these cycles are strongly dependent on the availability of water and energy. These interactions shape future terrestrial biosphere responses to global change. Here, we present a new terrestrial ecosystem model, QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system), which has been designed from scratch to allow for a seamless integration of the fully coupled carbon, nitrogen, and phosphorus cycles with each other and also with processes affecting the energy and water balances in terrestrial ecosystems. This new model includes (i) a representation of plant growth which separates source (e.g. photosynthesis) and sink (growth rate of individual tissues, constrained by temperature and the availability of water and nutrients) processes; (ii) the acclimation of many ecophysiological processes to meteorological conditions and/or nutrient availability; (iii) an explicit representation of vertical soil processes to separate litter and soil organic matter dynamics; (iv) a range of new diagnostics (leaf chlorophyll content; ¹³C, ¹⁴C, and ¹⁵N isotope tracers) to allow for a more in-depth model evaluation. In this paper, we present the model structure and provide an assessment of its performance against a range of observations from global-scale ecosystem monitoring networks. We demonstrate that QUINCY v1.0 is capable of simulating ecosystem dynamics across a wide climate gradient, as well as across different plant functional types. We further provide an assessment of the sensitivity of key model predictions to the model's parameterisation. This work lays the ground for future studies to test individual process hypotheses using the QUINCY v1.0 framework in the light of ecosystem manipulation observations, as well as global applications to investigate the large-scale consequences of nutrient-cycle interactions for projections of terrestrial biosphere dynamics.

1 Introduction

Past, present, and future changes in climatic conditions and atmospheric CO₂ concentrations affect terrestrial vegetation and soils (Hou et al., 2018; De Kauwe et al., 2013; Swann et al., 2016), which in turn induce biogeophysical and biogeochemical feedbacks to the atmosphere (Bonan, 2008; Friedlingstein et al., 2014; Zaehle et al., 2010). To predict the likely trajectories of terrestrial ecosystems under climate change and their climate feedbacks, it is important to develop and test advanced modelling tools for the terrestrial biosphere (Sitch et al., 2015). Global terrestrial biosphere models (TBMs) have evolved during the last decades alongside our understanding of soil and vegetation functioning (Bonan and Doney, 2018). Early vegetation models included only processes related to the carbon (C) and water cycle (e.g.  Sitch et al., 2003; Krinner et al., 2005; Dickinson et al., 2006; Raddatz et al., 2007; Clark et al., 2011). The recognition that nitrogen (N) has a pivotal role in the future dynamics and C storage of terrestrial ecosystems (Vitousek and Howarth, 1991; Hungate et al., 2003) has led to the development of a new generation of models that include a comprehensive representation of the N cycle (Thornton et al., 2007; Zaehle and Friend, 2010; Gerber et al., 2010; Wårlind et al., 2014). More recently, phosphorus (P) cycle processes have also been included in TBMs (Wang et al., 2010; Goll et al., 2017; Yang et al., 2014), as these can be important to understanding ecosystem response in naturally P-limited tropical regions (Yang et al., 2014) or in regions with large atmospheric N deposition and ensuing increases in P limitation (Peñuelas et al., 2013; Vitousek et al., 2010).

Although the number of element cycles and thereby processes considered in TBMs has increased in an effort to provide more realistic models (Zaehle and Dalmonech, 2011; Prentice et al., 2015), there is little consensus on how to represent many of these processes in a realistic but robust and computationally efficient manner. Often, small-scale soil and vegetation processes can lead to larger-scale feedbacks (e.g. Hararuk et al., 2015; Bradford et al., 2016) and therefore need to be included in sufficient detail in terrestrial biosphere models. Recent model–data synthesis studies with observations from Free-Air CO₂ Enrichment (FACE) experiments have revealed a number of fundamental issues in the way that state-of-the-art models represent plant growth (De Kauwe et al., 2014), nitrogen dynamics (Zaehle et al., 2014), and water responses (De Kauwe et al., 2013) to elevated CO₂. This highlights the need for a modular framework that allows testing multiple hypothesis for uncertain processes, thereby gaining an understanding on how process information and importance propagates from site to regional and global scales (Knauer et al., 2016; Peng et al., 2014; Thum et al., 2017).

One important obstacle to such an approach is the gradual development of terrestrial biosphere models, which implies that new features, such as processes describing the dynamics of the N or P cycle, have been added to existing carbon-cycle-only land surface model formulations. This evolutionary approach can result in a situation where assumptions that were made in earlier versions of the model are incompatible with the new assumptions, or that the old model structure cannot appropriately accommodate new structures, therefore limiting the ability to take new ecophysiological understanding into account. To overcome such issues, we have developed a new terrestrial biosphere model, QUINCY (QUantifying Interactions between terrestrial Nutrient CYcles and the climate system; Fig. 1), benefiting from advances in the understanding of global ecology. This new model is formulated based on the following:

• The recognition that plant growth may be limited by source (e.g. photosynthesis) or sink (growth rate of individual tissues, constrained by nutrients, temperature, and water availability) processes (Fatichi et al., 2013): we account for this decoupling by introducing a short-term (“labile”) and a long-term (“reserve”) storage pool for carbon and nutrients, which allow plants to adjust the carbon partitioning between growth, respiration, and storage according to nutrient availability, moisture stress, and temperature.

• The consideration of gradual changes in nutrient availability over time and their effects on the surface carbon, water, and energy exchanges: the model represents all biogeochemical and biogeophysical processes in the model at a half-hourly time step. Many processes, e.g. the response of photosynthesis, tissue stoichiometry, and fine root growth to nutrient shortage, are assumed to have a process-specific “memory” timescale, causing a lagged response to instantaneous variations in the environmental conditions. This includes a representation of the acclimation of both photosynthesis and maintenance respiration to the prevailing growth temperature (June et al., 2004; Atkin et al., 2014; Mercado et al., 2018). Together with a dynamic representation of photosynthetic capacity based on soil nutrient availability, this feature reduces the need for regionally defined plant functional types (e.g. boreal versus temperate type) needed to describe the vegetation in different regions.

• The recognition that biogeochemical processes in the soil (e.g. soil organic matter decomposition, nitrification, denitrification, and weathering) vary strongly within the soil profile (Ahrens et al., 2015; Koven et al., 2013): therefore, soil physics, soil biogeochemical processes, and their interaction with vegetation processes through plant nutrient/water uptake as well as root turnover are explicitly vertically resolved.

Figure 1 Schematic representation of the model structure. Ellipses: biogeochemical pools and other state variables; rectangles: biogeochemical processes; tetraethers: model input; solid green lines: carbon fluxes; solid dark red lines: nitrogen and phosphorus fluxes; solid black lines: carbon, nitrogen, and phosphorus fluxes; dotted blue lines: effects.

The QUINCY model contains entirely newly written code, although certain process representations are adapted from literature and previous models, including but not limited to OCN (Zaehle and Friend, 2010) and JSBACH3 (Roeckner et al., 2003). This new code approach allows for an appropriate separation of model infrastructure (e.g. memory allocation) and scientific code, and a better integration of model components. This permits us to include an internally consistent representation of ecosystem processes and in particular nutrient effects on plant photosynthesis, growth, and soil organic matter turnover through the inclusion of a common set of underlying hypotheses.

The aim of this paper is to present the basic structure and main features of the baseline version of this new model. We provide an assessment of the model's performance against a range of observations from global-scale ecosystem monitoring networks to demonstrate that the framework is capable of simulating ecosystem C fluxes in agreement with these observations given the simulated N, P, and water availability across a large gradient in climate and soil conditions, as well as across different plant functional types. To aid this understanding, we provide an assessment of the model's sensitivity to its parameterisation.

The information from this paper provides the background for future process-specific studies to test and improve process representation. In particular, we have formulated this model in a modular structure to facilitate the testing of multiple hypotheses for one particular process and to reduce the scope of the model (e.g. separately simulating canopy dynamics, vegetation dynamics, or soil biogeochemistry) to study the impact of particular processes (see the Supplement for details). Additionally, we have added a number of important diagnostics that will allow further ways of model evaluation. For instance, the leaf chlorophyll content is explicitly modelled from the leaf N to improve the simulated light response of photosynthesis throughout the canopy (Kull and Kruijt, 1998), but this will also allow for novel ways to compare the simulation results with in situ and remotely sensed observations in the future. The model also keeps track of the isotopic composition of C (in terms of ¹³C and ¹⁴C) and N (¹⁵N), which will allow to make use of natural abundance data and isotope tracer studies in the future.

2 Methods

2.1 Model description

QUINCY model v1.0 (rev. 1996) represents the coupled biogeochemical cycles of carbon (C), nitrogen (N), and phosphorus (P) in terrestrial ecosystems, as well as their interaction with the terrestrial water and energy balance (Fig. 1). The model traces the flow of these elements as well as C and N isotopes (¹³C and ¹⁴C, and ¹⁵N) through vegetation and soil at a half-hourly time step. The model considers eight broadly defined plant functional types (PFTs; see Table 1), representing differences in leaf type (herbaceous, broad-leaved, coniferous), phenology (evergreen, rain green, and summer green), growth form (grasses and trees), and photosynthesis type (C₃ versus C₄ photosynthesis) with a set of traits describing time-invariant functional differences across types (see Table S7 in the Supplement). Different from other TBMs, certain “soft” plant functional type characteristics, such as the photosynthetic capacity of leaves or the temperature sensitivity of photosynthesis and respiration acclimate to local environmental conditions, lead to a smoother transition of ecosystem functioning across ecoclimatic gradients. The number of PFTs is not embedded into the code structure and therefore new PFTs can be easily added if sufficient data to parameterise these are available.

Table 1 Description of plant functional types used in the model.

A brief overview of the model structure and key processes is given in this section. The detailed mathematical description of the model processes can be found in the Supplement. Although conceived to be coupled to the land surface scheme of an Earth system model, the model is currently applied stand-alone at the ecosystem scale. The stand-alone version requires half-hourly short- and longwave radiation, air temperature, precipitation and snowfall, wind velocity, atmospheric CO₂, ¹³CO₂, ¹⁴CO₂ mole fractions, as well as rates of NH_x, NO_y, and PO₄ deposition as time-dependent driving data. In addition, the model requires input describing the geographical coordinates, vegetation type, and soil physical and chemical parameters (texture, bulk density, rooting and soil depth, as well as inorganic soil P content).

2.1.1 Vegetation processes

Expanding on the concept employed by the Lund–Potsdam–Jena (LPJ) model (Sitch et al., 2003) and OCN (Zaehle and Friend, 2010), vegetation is characterised by an average individual representative of a PFT, composed of three fast-lived structural tissue types (pools: leaves, fine roots, and fruits), a respiring non-structural pool (labile), as well as a seasonal, non-respiring, and non-structural storage pool (reserve). In the case of tree vegetation types, three longer-lived structural tissue types (sapwood, heartwood, and coarse roots) are additionally represented, as are stand characteristics (height, diameter, and individuum density).

Building on Zaehle and Friend (2010), photosynthesis is calculated for sunlit and shaded leaves separately, explicitly taking the vertical canopy gradient of light, foliar chlorophyll, and photosynthetic N into account (Kull and Kruijt, 1998), with extensions for C₄ photosynthesis (Friend et al., 2009) and the temperature response as in Bernacchi et al. (2001, Supplement Sect. S2). Different from OCN, temperature acclimation of photosynthesis is modelled to acclimate to air temperature according to Friend (2010). Photosynthetic parameters, including chlorophyll content, further are assumed to depend on leaf N concentration (Friend et al., 1997). As foliar mass and N concentrations both respond to soil N availability (Vicca et al., 2012; Hyvönen et al., 2007; Meyerholt and Zaehle, 2015), this causes canopy photosynthesis to be directly affected by soil N availability in addition to meteorological variables. The stomatal conductance is coupled to the calculation of photosynthesis according to the Ball–Berry formulation (Ball et al., 1987; Knauer et al., 2015). Different from Zaehle and Friend (2010), QUINCY v1.0 accounts for stomatal and biochemical limitations by soil moisture (Egea et al., 2011, Supplement Eq. S18). In addition, photosynthesis can become downregulated due to sink limitation, when nutrient or water availability or low temperatures limit growth and cause accumulation of photosynthates in the labile pool (Fatichi et al., 2013; Hartmann et al., 2018, see Fig. S1 in the Supplement for an example). Plant nutrient uptake is assumed linearly dependent on fine root biomass density for each soil layer and follows a Michaelis–Menten parameterisation to simulate the effects of soil soluble NH₄, NO₃, and PO₄ concentrations (Zaehle and Friend, 2010, Supplement Sect. S4.5). As in Zaehle and Friend (2010), plant internal nutrient demand can up- or downregulate uptake, but rather than relying on instantaneous plant demand, the response of plant nutrient uptake to plant demand is modelled as a lagged response (of a few days) to balance short-term fluctuations in photosynthesis and soil nutrient availability and to represent memory effects in the plant's control of its nutrient uptake (Fig. S2).

As in LPJ and OCN, maintenance respiration is calculated for each tissue type as a function of its N content and tissue temperature (Sect. S3.3, Lloyd and Taylor, 1994). However, a novel feature of QUINCY v1.0 is that it does account for the observed acclimation of basal respiration rates to growth temperature (Atkin et al., 2014). A second new feature is that the costs for actual nutrient transformation, which are dependent on the rate of uptake and source of nitrogen (Eq. S24, Zerihun et al., 1998), are included in the root maintenance respiration term. While maintenance respiration takes priority over growth, under severe C deficit after the downregulation of photosynthesis, the maintenance respiration and nutrient uptake respiration can become downregulated by the meristem activity.

Different from earlier model approaches that included a labile pool (Zaehle and Friend, 2010), plant growth is modelled as the balance of source processes (photosynthesis, nitrogen and phosphorus uptake) and the capacity of the plant to create sinks (production of new biomass tissue, respiration, and storage; Sect. S3.2). The capacity of the plant to generate sinks is dependent on (i) the activity of the plant's meristem controlling the usage of the labile pool for growth, which is reduced at low air temperature and/or soil moisture (Fatichi et al., 2013); (ii) the partitioning of its labile resource to new tissue growth and the availability of nitrogen and phosphorus to create these tissues (Zaehle and Friend, 2010); (iii) the need for longer-term storage production (Fisher et al., 2010); and (iv) respiration for maintenance, which is given priority over growth of new tissue. Short-term fluctuations in the balance between carbon acquisition through photosynthesis and nutrient acquisition by roots are balanced by the timescale of the labile pool turnover. Seasonal and longer-term fluctuations are balanced to some extent by the build-up and use of the reserve pool. The flow of compounds between the labile and reserve pools (Sect. S3.6) is described by a set of pull (from the reserve pool) and push (to the reserve pool) mechanisms attempting to simulate the signalling related to the beginning of the growing season (phenological pull), the need for maintenance and growth (maintenance pull), as well as to provide resources for the next growing season or resilience against interannual variability (reserve push).

Tissue growth follows a set of allometric equations (Shinozaki et al., 1964; Sitch et al., 2003; Zaehle and Friend, 2010), describing size-dependent relationships of foliar, fine root, coarse root, sapwood, and fruit biomass (Sect. S3.4). The allometric relationship between leaves and fine roots responds to N, P, and water limitation and therefore increases uptake capacity under nutrient limitation. The C:N and N:P ratios of the leaves and fine roots are flexible within empirical bounds and respond dynamically to an imbalance between C and N, or N and P availability, respectively (Sect. S3.5). The allometric equations and the N and P content of each tissue give the stoichiometric constraint of plant growth. Tissue turnover is considered at two timescales: (i) the rapid breakdown of enzymes (N) and associated RNA (P) from tissue, and their ensuing replacement from the labile N and P pools (recycling; at the timescale of days; Zaehle and Friend, 2010), which allows for seasonal changes in tissue nutrient concentrations; and (ii) the senescence of entire tissues and their conversion to litter, including the resorption of nutrients from foliage (Sect. S3.9). The tissue senescence is calculated given a fixed turnover time for each tissue, with the exception of the leaves, where turnover is determined by the PFT-specific phenological timing (Sect. S3.8).

Stand-level vegetation dynamics are represented through the establishment and mortality of an average representative individual of the entire tree populations, following largely Sitch et al. (2003) and Zaehle and Friend (2010, see Sect. S3.10). Mortality is simulated as density dependent, constraining the number and size of individuals, as well as a dynamic background mortality related to the overall growth efficiency of the plant. (Re-)establishment and recruitment is dependent on meteorological and space constraints, but a new aspect of the re-establishment representation is that the model explicitly takes account of the available seed pool at the forest/grassland floor, which depends on the vegetation's ability to grow and produce fruit.

2.1.2 Soil biogeochemical processes

The turnover and formation of soil organic matter (SOM) follows, with some modifications, the CENTURY approach (Parton et al., 1993, Sect. S4): the turnover of litter and SOM pools are all calculated using first-order kinetics with temperature and moisture dependence. The litter stoichiometry is determined by the stoichiometry of senescent plant tissue and tissue-specific allocation fractions. The stoichiometry of the fast SOM pool adjusts to available inorganic nutrients, whereas the slow SOM pool is assumed to have fixed $\mathrm{C}:\mathrm{N}:\mathrm{P}$ stoichiometry. As a result, the decomposition of litter may become N and P limited, leading to the accumulation of an organic litter layer and reducing the rate of N and P cycling in the ecosystem.

The soluble inorganic NH₄, NO₃, and PO₄ are assumed to be available for plant uptake and microbial immobilisation, as concluded by White (2006) and many others. Plant uptake and SOM decomposition compete for these nutrients based on their respective demand and uptake capacity, which vary in time and with soil depth. In the aerobic part of the soil, NH₄ is oxidised to NO₃ through nitrification and in the anaerobic part of the soil NO₃ is reduced to N₂ through denitrification (Zaehle et al., 2011). Both processes produce NO_y and N₂O as byproducts (Sect. S4.7). The representation of soil inorganic P dynamics (Sect. S4.8) follows the concept from Yang et al. (2014) and Wang et al. (2010). The available soluble PO₄ is exchanged via ad-/desorption with the mineral surface, leached with soil water and replenished by atmospheric deposition, weathering, and biomineralisation of P that is bound in the slow SOM pool. The adsorbed P is further absorbed into the soil matrix and gradually becomes an inactive form of P (occluded P;  Walker and Syers, 1976).

Different from the original CENTURY-style models, the vertical soil profile of biogeochemical pools and processes is explicitly represented, with exponentially increasing layer thickness with increasing soil depth. For each of these layers, the model represents different organic pools (metabolic, structural, and woody litter, as well as fast and slow overturning soil organic matter (SOM)), as well as inorganic pools of N and P (soluble inorganic NH₄, NO₃, as well as NO_y, N₂O, and N₂, soluble inorganic PO₄, adsorbed PO₄, occluded PO₄, and primary PO₄). For each soil layer, temperature and moisture are calculated based on soil physical characteristics, and the transport and atmospheric exchange of energy and water (Sect. S6.3). As an example of the benefit of QUINCY's modular approach, Fig. S3 shows that the explicit representation of the vertical soil profile, compared to a zero-dimensional, lumped soil approach, has little effect on the seasonal course of heterotrophic respiration. However, it does affect the simulated nutrient dynamics because of the explicit separation of a nutrient immobilisation in the litter dominated layers from the gross-mineralisation-dominated soil layers with a proportionally higher content of soil organic matter.

2.2 Data for model evaluation

We evaluate simulated diurnal and seasonal patterns of gross primary production and net ecosystem carbon exchange (GPP and NEE, respectively) at a subset of FLUXNET sites (see Table A1) available from the FLUXNET La Thuile synthesis dataset (NOA, 2007). The half-hourly CO₂ fluxes have been measured and processed using standard procedures (Papale et al., 2006) as well as gap-filling and partitioning algorithms (Reichstein et al., 2005). We further evaluate the mean annual estimates of gross and net primary production (GPP and NPP, respectively) and their ratio, the carbon-use efficiency (CUE = NPP ∕ GPP) from a global forest database (GFDB; Vicca et al., 2012; Campioli et al., 2015). We further evaluate the simulated growing-season plant water-use efficiency, i.e. the ratio of plant water loss to carbon uptake, by comparing a proxy, foliar isotope discrimination of ¹³C (Medlyn et al., 2017) against data reported by Cornwell et al. (2018a, b). Figure 2 provides an overview of the sites used for evaluation.

Figure 2 Locations of the sites used for model evaluation. Black stars: FLUXNET sites; circles: GFDB sites with circles having colours corresponding to different PFTs (see Table 1).

2.3 Model setup

2.3.1 Boundary conditions and meteorological forcing

The QUINCY model is applied at site scale for the dominant PFT at each site. For the simulations at the FLUXNET sites, the site-specific observed meteorological data were used (NOA, 2007). For the GFDB sites, where site-level meteorology is not readily available, daily meteorological data for 1901 to 2015 were taken from Climate Research Unit – National Centers for Environmental Prediction (CRUNCEP), version 7 (Viovy, 2016) and disaggregated to the model time step using the statistical weather generator as in Zaehle and Friend (2010). The annually changing CO₂ mole fraction was obtained from Le Quéré et al. (2018), and the time series of ¹³C and ¹⁴C were obtained from Graven et al. (2017). The time series of N deposition for each site was taken from Lamarque et al. (2010) and Lamarque et al. (2011). For the P deposition model, estimates of nutrient and dust fluxes from Brahney et al. (2015) and Chien et al. (2016) were used.

Soil physical properties (volumetric heat capacity and conductivity, water content at saturation, field capacity, and wilting points, as well as parameters describing the soil water retention curve) are derived from soil texture according to Saxton and Rawls (2006). Where available, texture data are taken from site observations. Alternatively, they are obtained from the nearest grid cell of the SoilGrids dataset (Hengl et al., 2017). The rooting depth has been obtained from Jackson et al. (1996).

The vegetation is initialised as bare ground with a small amount of seed-bed mass to start vegetation growth. The soil organic matter profile is initialised with a default SOM content for each pool, 60 % of which is allocated to the uppermost layer, and the remaining 40 % is allocated to other soil layers in proportion to the assumed rooting profile. Previous tests have shown that the model is not very sensitive to the choice of the SOM initial conditions. The soil inorganic P pools of the first 50 cm are initialised using the soil P dataset by Yang et al. (2013), extrapolated to the whole soil profile assuming (i) a constant total soil inorganic P density (mol P kg⁻¹ soil) and (ii) an increasing fraction of primary P and decreasing fraction of labile P with increasing soil depth.

2.3.2 Spinup and model protocol

The QUINCY model is spun up with repeated meteorological forcing for each site to obtain an near-equilibrium soil and vegetation state. For the GFDB sites, meteorology and other atmospheric forcing (atmospheric CO₂, its isotopes, as well as N and P deposition) are taken by randomly drawing years from the period of 1901 to 1930. The same procedure is applied for FLUXNET sites, with the exception that meteorological forcing is only available for more recent years, and data from those years are used in random sequence. The P cycle is activated during the model spinup, but the more stable soil inorganic P pools, i.e. the primary, sorbed, occluded pools, are kept constant to ensure that the simulation initialises at the P status of Yang et al. (2013)'s dataset.

After spinup, the model is run for the period of 1901–2015 using the annual values for atmospheric CO₂, ¹³CO₂, and ¹⁴CO₂, as well as N and P deposition, and the meteorology of the respective year. For the comparison to GFDB data, forest stands are clear cut in the year provided by the database. Upon harvest, vegetation biomass is transferred to the litter pools, with the exception of the wood pool, of which a given fraction, $f_{\mathrm{harvest}}^{\mathrm{wood}}$, is removed from the site. In the case of FLUXNET sites, we retain the model output for the years covered by the available data at each individual site (Table A1). For the GFDB sites, we average over the period of 1995–2004 to remove the effects of interannual climate variability and to have modelled forest age close to the actual forest age of each site, as the precise year of measurement is not recorded in the database for all sites.

To test the effect of the simulated N and P availability, we applied the model for the FLUXNET sites with three different setups. Next to the full CNP version of the model described above (referred to hereafter as “CNP”), we performed a simulation (“CN”) in which the P concentration of the soil was kept at concentrations not limiting plant uptake or SOM decomposition, effectively removing the impact of the inorganic P sorption and weathering dynamics, and maintaining N:P stoichiometry at default values. We also added a version (“C”), in which, in addition, biological N fixation in soil (asymbiotic fixation) and in vegetation (symbiotic fixation) was calculated so as to avoid any N limitation of soil turnover and vegetation growth while keeping the C:N stoichiometry at default values.

2.4 Model sensitivity to its parameterisation

We further test the sensitivity of the model to its parameterisation using a hierarchical Latin hypercube sampling design (LHS; Saltelli et al., 2000; Zaehle et al., 2005) to assess the stability of the model with respect to changing parameter values. As many parameters have unknown value ranges and distributions, we simply vary each parameter between 90 % and 110 % of its default value given in the Supplement (Tables S1–S9), drawn with LHS from a uniform distribution. We first generate LHS samples for each module (corresponding to a section in the Supplement; results not shown) and evaluate the model output from these simulations in terms of long-term mean modelled GPP, leaf $\mathrm{C}:\mathrm{N}:\mathrm{P}$, net N and P mineralisation, as well as total vegetation and ecosystem carbon storage. Based on these simulations we select parameters (n=45), which have a strong effect on the model outcome to form a new set of 1000 LHS samples, which are used for the model sensitivity analysis presented in this paper. We measure parameter importance as the rank-transformed partial correlation coefficient (RPCC) to take account of potential non-linearities in the relationship between parameters and model output (Saltelli et al., 2000; Zaehle et al., 2005).

3 Results

3.1 Simulated diurnal to seasonal timescale carbon dioxide fluxes at selected FLUXNET sites

We first compare the simulated mean monthly diurnal as well as the mean seasonal cycle of GPP at four forest sites, representing the major tree plant functional types, with in situ observations (Fig. 3; see Fig. S4 for a comparable evaluation of the latent heat flux). The model simulations agree reasonably well with the observations in that diurnal and seasonal timing and overall magnitude of fluxes are generally consistent. Importantly, while the N and P cycles dampen the overall magnitude of the fluxes at the beech forest site (temperate broad-leaved deciduous; FR-Hes) and the eucalyptus forest site (temperate broad-leaved evergreen; AU-Tum), they do not affect the shape of either the diurnal or seasonal cycles. In other words, the effect of diurnal and seasonal variability in soil nutrient availability is buffered through the labile and reserve storage pools in the vegetation, such that it affects vegetation gross carbon uptake only via slow processes such as foliar nutrient and allocation changes but has no effect on variability at the daily to weekly timescale. This is also demonstrated in the LAI values that are influenced by the long-term dynamics. For instance, at FR-Hes and AU-Tum the LAI is lower at the nitrogen- and phosphorus-dynamics-enabled version than with the C-only version (6.0 compared to 6.2 m² m⁻² at FR-Hes, 3.7 compared to 5.9 m² m⁻² at AU-Tum), whereas it does not have a notable effect for the needle-leaved evergreen site of FI-Hyy.

Figure 3 Simulated and observed mean monthly diurnal (a, b, c, d) and seasonal (e, f, g, h) cycles of GPP at four FLUXNET sites (FI-Hyy, FR-Hes, AU-Tum, and BR-Ma2; see Table A1) representing the major QUINCY PFTs (NE, BS, TeBE, and TrBE, respectively; see Table 1). “Obs” corresponds to micrometeorological observations. “C”, “CN”, and “CNP” refer to the model simulations with C, C&N, and C&N&P options enabled. Seasonal cycles have been smoothed by a 16 d running mean.

There is little overall effect of the nutrient cycles on the correlation between half-hourly simulated and observed GPP at FI-Hyy (all versions: r²=0.76), FR-Hes (all versions: r²=0.63), and BR-Ma2 (all versions: r²=0.67), and only a minor reduction for AU-Tum (r² 0.70 versus 0.65 for C and CN(P), respectively). The largest decrease of the model's root mean square error (RMSE) took place at FR-Hes, where introducing N and P cycling reduced the RMSE by nearly 20 % to 5.17 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. For FI-Hyy, the RMSE and the improvement were 2.42 and 0.11 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. For AU-Tum, the RMSE and worsening of RMSE were 5.44 and 0.36 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. At BR-Ma2, the RMSE and the improvement were 7.78 and 0.09 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. At the tropical evergreen forest site BR-Ma2 (Fig. 3d, h), the different versions performed similarly. The lower GPP in the daily observation-based estimate between days of year 150 and 200, which are not replicated by the model, is largely caused by lower observation-based GPP at the start of the time series during this time of year in 2003. There was no obvious cause for this behaviour visible in the meteorological drivers. The latent heat fluxes were generally well simulated at these four sites, with the exception of large summertime overestimation (73 %) at FR-Hes (Fig. S4).

Table 2 summarises the observed and simulated GPP, total ecosystem respiration (TER), and foliar C:N with the CNP version of the model for these four FLUXNET sites, calculated for the time period of the flux observations (Table A1). The annual GPP from the CO₂ flux observations is in line with the modelled values for FI-Hyy and FR-Hes, while productivity at BR-Ma2 is somewhat overestimated and strongly underestimated at AU-Tum. The TER is somewhat overestimated for FI-Hyy and FR-Hes, and strong overestimation occurs at BR-Ma2. However, the unusually large discrepancy between observed GPP and TER fluxes at BR-Ma2 suggests that these estimates should be considered with caution. At AU-Tum, the TER is underestimated. Simulated foliar C:N was within the observational range reported in the La Thuile database (NOA, 2007) for all four sites. Observations for nutrient availability and fluxes are not consistently available across these sites.

Table 2 Simulated and observed GPP (in $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), total ecosystem respiration (TER in $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$), and leaf C:N (gC gN⁻¹) at four FLUXNET sites, as in Fig. 3. Values are averaged over the years where CO₂ flux observations are available.

The Taylor plots of the half-hourly simulated and observation-based time series of GPP and NEE (Fig. 4) demonstrate that the overall model performance noted in Fig. 3 holds across a wider set of sites, spanning a larger climatic gradient and also including grassland sites. The Taylor plot reports results obtained with the CNP version of the model; the plots for the C and CN version were very similar in terms of their correlation and moderately different in terms of the ratio of standard deviations. The simulated GPP shows slightly better performance over the different sites than NEE. Considering that (i) the model has not been parameterised specifically for each site and (ii) eddy-covariance observations at this timescale are subject to considerable fluctuations, the model performance is acceptable in terms of the correlation. While there is a wide spread in the ratio of the simulated to observation-based standard deviation, mostly associated with site-based differences in the seasonal maximum of the flux, there is no systematic model bias for any PFT.

Figure 4 Taylor plots for GPP and net ecosystem exchange (NEE) showing the model–data agreement for the FLUXNET sites used in this study (Table A1), separated according to the dominant plant functional type (Table 1). The standard deviation was normalised against the standard deviation of the observations for the corresponding variable. The grey lines correspond to the Euclidean distance from the point of perfect model–data agreement, where both the normalised standard deviation and the correlation coefficient are 1.0, shown as a black star in the figure.

The inclusion of the N cycle (difference between C and CN model versions) reduces the RMSE between simulated and observation-based estimates of GPP for all PFTs apart from the tropical broad-leaved evergreen PFT (Fig. 5). The largest reductions of the RMSE between the different versions occur for TeBE (34.0 %) and TeH (41.6 %). There is no strong effect of including the P cycle on simulated GPP, and in consequence there is no difference in the RMSE values of the CN and CNP versions, even for the tropical broad-leaved evergreen forest sites, which are in regions typically associated with P limitation of biomass production.

Figure 5 Root mean square error (RMSE) for simulated GPP by dominant plant functional type (see Table 1). The number of sites used in the calculation of the PFT-specific RMSE value (n) is shown above the bars for each PFT. The error bars denote the standard deviation of the RMSE values of the different sites within each PFT. “C”, “CN”, and “CNP” refer to the model simulations with C, C&N, and C&N&P options enabled.

3.2 Cross-site analysis of long-term mean GPP and NPP

Simulated GPP across all the GFDB sites (n=434), covering a range of tree PFTs, agrees reasonably well with data (r²=0.55, RMSE =560 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) (Fig. 6a), despite an underestimation of GPP for some sites, both in tropical forest ecosystems and needle-leaved evergreen forests. The data shown are restricted to forests older than 20 years to avoid ambiguities in the model–data comparison as the exact method of forest regeneration after disturbance is not known and the model does not describe stand-level dynamics of very early succession forests that often have very high observed productivity. The model is also able to replicate the observations from low to modest NPP values albeit with a higher relative RMSE (Fig. 6b) (r²=0.35, RMSE =305 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). While there is no systematic over- or underestimation of NPP with PFT or magnitude of NPP, it is noteworthy that the highest NPP value simulated for this dataset was only 1441 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, whereas for a few sites, observation-based estimates reached nearly 2000 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. The model correctly predicts the range of CUE, i.e. the ratio of net to gross primary production, where the observed CUE values ranged between 0.21 and 0.76 g g⁻¹, while the simulated values were between 0.19 and 0.61 g g⁻¹. The model also predicts correctly a lower CUE in tropical compared to temperate/boreal ecosystems. However, the model is not able to explain the within-PFT variability of the observed CUE values in the dataset (Fig. 6c).

Figure 6 Simulated versus observed GPP (a), NPP (b), and carbon use efficiency (CUE = NPP ∕ GPP) (c) at GFDB sites against observations. For PFT abbreviations, see Table 1. Note that observations of GPP and NPP are not consistently available for all sites.

Underlying these results are the emergent sensitivities of ecosystem processes to the model drivers. Figures S5 and S6 show that the large-scale gradients of simulated GPP and NPP are largely controlled by mean annual temperature (MAT), with an additional clear effect of mean annual precipitation. For temperate and boreal forests, also N deposition is positively related to GPP and NPP. While the acclimation of photosynthesis to growth temperature does matter at the diurnal timescale particularly on cloudy days (Fig. S7), the positive and negative effects cancel each other out at the long-term annual mean and therefore play no role in the simulated GPP–MAT relationship. Despite a notable decline of maintenance respiration in higher temperatures due to the acclimation of respiration to growth temperature (see Fig. S8), the spatial trend in simulated CUE is predominantly driven by MAT (Fig. S9). Trends apparent in the CUE response to MAP and N deposition are small and likely confounded by the MAT response. However, the general tendency of PFT-specific foliar C:N to decline with increasing N deposition (Fig. S10), something that is expected to happen based on observational evidence (Hyvönen et al., 2007; Meyerholt and Zaehle, 2015), suggests that the linearly increasing maintenance respiration rates per unit tissue drives whole-plant maintenance respiration up stronger than the saturating effect of increased leaf N on GPP, therefore reducing CUE with increasing N availability, contrary to expectations (Vicca et al., 2012). The root-to-shoot ratio did not show strong dependency to any of these environmental variables (data not shown), probably contributing to the lack of a decline in CUE with declining N availability. Generally, PFT-wise foliar stoichiometry appears to be more strongly influenced by N deposition than MAT or MAP. This result occurs despite a clear dependency of the simulated biological N fixation on temperature and precipitation (Fig. S11), where the simulated range and response to climate drivers are consistent with the available evidence suggesting the highest N fixation in hot and moist environments (Cleveland et al., 1999, 2013). In addition, in agreement with recent evidence, enhanced atmospheric N input suppresses simulated N fixation (Zheng et al., 2019).

The comparison of leaf Δ¹³C based on observations (Cornwell et al., 2018b) and the GFDB simulations (Fig. 7) provides another look into the underlying processes of GPP evaluation, as they are documenting the simulated and observed trends in water-use efficiency with water availability. Both observations and model showed lower Δ¹³C in drier conditions, indicating an effect of seasonal water availability in overall drought stress on photosynthesis. Comparing PFT-wise mean values, constrained to estimates of ET∕PET >0.25 to remove the effect of modelled outliers under high water limitation, the average for broad-leaved deciduous forest (BS) was 20.88 ‰ for the observations and 20.67 ‰ for the modelled values. For the needle-leaved deciduous (NS) PFT, the average values were very similar (observations 20.29 ‰ and simulations 20.12 ‰). Somewhat larger discrepancies occurred for the needle-leaved evergreen PFT, where the average of the observations was 20.00 ‰ and the average of the simulation results was 19.59 ‰. For the broad-leaved evergreen forests, there was a large deviation between the observed (22.68 ‰) and modelled values (20.55 ‰), but also the sample sizes were very different (n=1617 for observations and n=54 for simulations). For the modelled values, the difference between needle-leaved deciduous and evergreen forests was more pronounced than in the observations, but the sample size of modelled NS sites was very small (n=6).

Figure 7 Observed and simulated leaf Δ¹³C. Observations are taken from the global database of Cornwell et al. (2018b); simulations are derived from the GFDB set of sites (see Fig. 2). As there is no direct correspondence between data and model location, the data are plotted as a function of climatic water deficit, measured as the ratio of actual to potential evapotranspiration (ET∕PET). For model PFT abbreviations, see Table 1.

3.3 Model sensitivity analysis

The parameter sensitivity study reveals that from the many parameters introduced in this model (see the Supplement), the model output is sensitive only to a limited amount of these parameters. The role of specific parameters varies across sites (Table 3) and with the output quantity of interest (Table S1). A number of interesting things can be noted when looking at which parameters dominate across sites (Table 3). In line with reports from other studies for other models (e.g. Zaehle et al., 2005; Zaehle and Friend, 2010; Wramneby et al., 2008; Massoud et al., 2019), the first set of parameters is related to those affecting the photosynthetic nitrogen-use efficiency, which directly affect the C assimilation per unit leaf area. In QUINCY v1.0, this is mostly related to the magnitude of leaf N assumed to be structural N ($k_{\mathrm{0}}^{\mathrm{struc}}$), the scaling of leaf N to leaf chlorophyll content ($k_{\mathrm{fn}}^{\mathrm{chl}}$, $a_{\mathrm{chl}}^n$), and the efficiency of RuBisCO-limited carboxylation rate per unit N ($v_{\mathrm{cmax}}^n$).

Table 3 The 10 most important parameters (P) determining model uncertainty at four FLUXNET sites (FI-Hyy, FR-Hes, AU-Tum, and BR-Ma2, as in Figs. 3 and 8) and reference to the respective parameter description table in the Supplement (T). Overall importance of parameters is measured by first calculating the RPCC for each output variable and then calculating the mean of the absolute RPCC values across GPP, net N∕P mineralisation, vegetation, and ecosystem C, as well as leaf $\mathrm{C}:\mathrm{N}:\mathrm{P}$, weighted by the uncertainty contribution of these model outputs. The parameter names are denoted by asterisks: ^* refers to photosynthesis-related parameters; ${}^{**}$ to soil-biogeochemistry-related parameters; ${}^{***}$ to parameters related to vegetation growth and dynamics. See also Table S9 for a detailed breakdown of RPCCs according to output variable.

The second set of parameters to which the model is sensitive, unlike other carbon-only models, relates to parameters determining the rate of litter and soil organic matter turnover (T_opt,decomp, ${\mathit{\tau }}_{\mathrm{slow}}^{\mathrm{base}}$), the efficiency of microbial SOM processing (${\mathit{\eta }}_{\mathrm{C},\mathrm{litter}\to \mathrm{fast}}$, ${\mathit{\eta }}_{\mathrm{C},\mathrm{fast}\to \mathrm{slow}}$), and the C:N of the slow decomposing soil organic matter (${\mathit{\chi }}_{{\mathrm{SOM}}_{\mathrm{slow}}^{X:Y}}$). All these parameters control the rate at which nitrogen and phosphorus are released by SOM decomposition, which in the QUINCY 1.0 model is important for the growth of leaves and associated photosynthesis. Such a dependency is to be expected by a nutrient-cycle-enabled model and this is clearly documented in the dependence of simulated GPP on simulated net N and P mineralisation (Fig. S12). It is important to note that for the effect of uncertainty in the temperature optimum of decomposition, T_opt,decomp, it is not the direct effect of the optimum point that matters, since it is larger than common soil temperatures for most of the sites, but a side effect of changing the default turnover time of litter and SOM along with the change in optimum temperature. Consistent with expectation, the rate of soil turnover is more important in cold that warm environments. Soil-based parameters are less important that photosynthetic parameters for the tropical site compared to the three non-tropical sites.

The third set of parameters identified as sensitive can be characterised as being related to vegetation growth and dynamics. On one hand, these are parameters dominating the carbon allocation within plants (k_latosa, sla), also as noted by Zaehle et al. (2005) and Massoud et al. (2019) for other models, as well as the tissue stoichiometry (${\mathit{\chi }}_{\mathrm{leaf}}^{X:Y}$, ${\mathit{\chi }}_{\mathrm{leaf}}^{X:Y}$). Importantly, also the scaling exponent k_rp, relating diameter to crown area, and thereby determining the space constraint for the establishment and mortality of tree PFTs, plays some role in some of the sites due to its effect on LAI and total vegetation biomass. At the tropical BR-Ma2 site, despite temperature acclimation, the basal maintenance respiration rate for leaves and fine roots $f_{\mathrm{resp},\mathrm{maint}}^{\text{non-woody}}$ is also a sensitive parameter, as is the scaling parameter for limiting nutrient uptake under high nutrient availability ($K_{\mathrm{demand}}^{\mathrm{half},X}$).

The parameter sensitivities investigated here propagate to uncertainty in the model output. The interquartile range of the model output (Fig. 8) is well constrained and centred around the results of the standard parameterisation reported in Table 2. Extreme parameter combinations cause larger variations but do not fundamentally change the model behaviour. Including wider parameter ranges for some more uncertain parameters will likely affect the absolute mean and uncertainty range of particular model output. Therefore, the results shown in Fig. 8 should not be interpreted as an estimate of parameter-based model uncertainty. However, the results do suggest that despite introducing a complex ecosystem model with many non-linear equations and parameters, the model does not reveal fundamental model instabilities.

Figure 8 Median, interquartile range, and absolute range of simulated GPP, net mineralisation of NH₄ (${\mathrm{\Phi }}_{{\mathrm{NH}}_{\mathrm{4}}}$) and PO₄ (${\mathrm{\Phi }}_{{\mathrm{PO}}_{\mathrm{4}}}$), leaf C:N and N:P, as well as vegetation C (veg. C) and total ecosystem C (up to a depth of 1 m) obtained for the four FLUXNET sites in Fig. 3 using Latin hypercube sampling (n=1000) for 45 parameters. Values have been normalised to the ensemble mean, given as the number for each site and variable to improve readability. The numbers below each individual box correspond the median value for each variable: GPP in units $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, ${\mathrm{\Phi }}_{{\mathrm{NH}}_{\mathrm{4}}}$ in units $\mathrm{gN}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, ${\mathrm{\Phi }}_{{\mathrm{PO}}_{\mathrm{4}}}$ in units $\mathrm{gP}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, and vegetation and total carbon in units kgC m⁻².

By varying parameters, we essentially generate an ensemble of simulations with identical climate and atmospheric nutrient input, but different soil nutrient-cycle rates and vegetation nutrient-use efficiencies. At the example of the broad-leaved deciduous forest site FR-Hes, the model shows, as expected, clear dependencies between the rates of net N and P mineralisation and GPP, which also propagate to the simulated carbon stock in vegetation (Fig. S12) underlying the uncertainty ranges reported in Fig. 8. For this example, there is also a clearly negative correlation between the net N mineralisation and foliar C:N, which gives rise to a negative covariation of GPP and leaf C:N. Interestingly, such a correlation does not exist for foliar N:P, probably due to a lack of effect of foliar P concentrations on photosynthesis.

4 Discussion

Land surface models with coupled carbon–nitrogen and phosphorus cycles have been published before (Goll et al., 2017; Wang et al., 2010; Yang et al., 2014). The QUINCY v1.0 model introduced in this paper distinguishes itself from these models in (i) that it seamlessly integrates nutrient dynamics with carbon, water, and energy calculations (e.g. there is no difference in time stepping between biogeophysical and biogeochemical processes) and (ii) that there is an explicit decoupling of the growth processes from C supply due to photosynthesis by accounting for temperature, moisture, and nutrient constraints in the formation of new tissue, as opposed to the simpler treatment of these dynamics in Zaehle and Friend (2010) and applied in Goll et al. (2017). The purpose of this paper is to provide a background for future studies on the effect of coupled biogeochemical cycles at large scales as well as on testing the effect alternative representations of important nutrient-cycle-related processes, for which the model is suitable due to its modular design. A detailed evaluation of all new submodules of the model is beyond the scope of this paper.

Overall, QUINCY v1.0 performs well in terms of diurnal, seasonal, and annual carbon fluxes across a range of ecosystem monitoring flux sites, spanning a large geographic gradient and a variety of ecosystem types. Despite the lack of any site-specific parameterisation, there is no systematic model bias for any single PFT for either the FLUXNET (Figs. 3 and 4) or the GFDB (Fig. 6) analyses. Including a coupling of the carbon cycle to representations of the nitrogen and phosphorus cycle affects long-term average productivity through its effects on photosynthetic capacity (changed foliar nitrogen concentration) and leaf area (as a result of the changes in root-to-shoot allocation). However, short-term diurnal or subseasonal shortage of soil nitrogen or phosphorus does not directly and instantaneously affect simulated productivity because of the buffering introduced through the labile and reserve storage.

Model predictions generally show an improvement in the prediction of GPP with the addition of the N cycle. This occurs despite a slight low bias in simulated foliar nitrogen concentration compared to the mean values of the TRY database (BS: 2.1 (0.6–3.2) % versus 1.7±0.3, TrBE∕TeBE: 1.7 (0.5–3.9) % versus 1.7±0.2, NE: 1.3 (0.5–1.9) % versus 1.0±0.4, NS: 1.9 (0.9–3.0) % versus 1.8±0.2 observed and simulated, respectively, Kattge et al., 2011). This comparison is indicative only, since there is a geographic mismatch between the spatial coverage and distribution of the modelled sites and the data in the TRY database. Note that these estimates are not fully independent because the minimum and maximum ranges of permissable foliar N were derived from Kattge et al. (2011).

The P-enabled version does not show additional improvement in simulating the regional differences in productivity. This is partly due to the fact that the sites with sufficient information on GPP and NPP available for model evaluation are as a majority located in temperate and boreal regions, where the main limiting nutrient is likely to be nitrogen (LeBauer and Treseder, 2008). However, even for the few tropical sites included in the analysis, where in general terms limitation of productivity by P would be expected (LeBauer and Treseder, 2008; Vitousek et al., 2010), we do not find any significant effect of the simulated P cycle. An important contributor is the fact that in QUINCY v1.0, the vertically explicit description of SOM dynamics permits very efficient recycling of organic P in the litter layer, where low mineral sorption capacities lead to efficient P uptake of soil microorganisms and vegetation. The efficient recycling in the litter layer is partially due to the biochemical mineralisation flux. Further observations are required in order to constrain this flux better (Reed et al., 2011). One further important factor is the unknown initial soil concentration of plant-available P, as well as uncertainties in the rate of P weathering. Improving the understanding of P availability across the globe remains an important challenge for the representation of P effects on productivity (Wang et al., 2010; Goll et al., 2017), requiring detailed observations of soil and vegetation P dynamics and manipulation experiments (e.g. Medlyn et al., 2016). In addition, similar to other models (Goll et al., 2017; Yang et al., 2014), QUINCY v1.0 currently assumes that P limitation solely affects productivity due to a stoichiometric constraint on growth, while other mechanisms may also play an important role in the acclimation of plant processes to different levels of P availability (Jiang et al., 2019).

The predicted GPP and NPP across a large range of climatic and biogeographic situations are in good agreement with observations and so is the average of the CUE (NPP ∕ GPP). However, the model does not reproduce the observed range in CUE in temperate and boreal forests, and the analysis of the main drivers of CUE shows that in the current model version, CUE is mostly a function of mean annual temperature. Vicca et al. (2012) suggested that this variance is associated with altered carbon partitioning and in particular increased belowground carbon allocation in response to nutrient shortage. The QUINCY v1.0 model simulates an increase of the root-to-shoot ratio with response to nutrient stress. However, the effect of this change on CUE is small, as increased root allocation not only decreases biomass production through increased allocation to higher turnover tissues, but also whole-plant-mass-specific respiration, given the implicit model assumption that mass-specific fine root respiration is smaller than leaf-mass-specific respiration. This inference is consistent with detailed observations at the FLUXNET site FI-Hyy (Ilvesniemi et al., 2009; Korhonen et al., 2013), where the QUINCY v1.0 model successfully simulated GPP and vegetation C storage (simulated/observed 7.0∕6.8 kgC m⁻², Ilvesniemi et al., 2009) but substantially overestimated the NPP (simulated/observed: 536∕242 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$). Additionally, the N uptake by the vegetation in addition to the N losses was in the same order of magnitude as the observations (Korhonen et al., 2013), suggesting that C partitioning rather than N availability is the source of the underestimation. Further causes of this model–data mismatch include alternative pathways of carbon partitioning not represented in the model (e.g. exudation), changes in tissue turnover times with changes in nutrient availability, and the magnitude of mass-specific autotrophic respiration. For instance, the estimated autotrophic respiration in FI-Hyy was 714 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Ilvesniemi et al., 2009), compared to 494 in the model $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, which may be the combined result of allocation, temperature response, and specific respiration rates, particularly below ground. The novel structure of the QUINCY v1.0 model that decouples photosynthesis from growth, and therefore permits to adjust carbon partitioning to different sinks at short timescales, will allow to investigate the effect of these processes in the future. Finally, the intra-PFT variation of observed CUE likely depends on other site-specific factors that affect site fertility, which are currently not adequately represented in the model, such as the effects of soil pH, site history, and species-level variability.

The comparison of simulated leaf ¹³C discrimination to observations recorded in the global dataset by (Cornwell et al., 2018b) suggests that the overall parameterisation of long-term mean leaf- and ecosystem-level water-use efficiency, derived from instantaneous leaf-level gas exchange measurements (Lin et al., 2015), broadly conforms with observations. Notably, the model also simulates the trend in discrimination with respect to changing water availability in accordance with observations, which suggests that the overall effect of moisture availability on water-use efficiency is appropriately considered by the model. Discrepancies between the observation-based estimates of water-use efficiency derived from leaf-level flux and isotope measurements have been noted before (Medlyn et al., 2017), and these may contribute to the remaining model–data mismatch. One possible reason for this mismatch may be the omission of mesophyll conductance in model formulations, which may induce systematic shifts in chloroplast [CO₂] and thereby affect the simulated value of leaf Δ¹³C without affecting leaf-level water-use efficiency (Knauer et al., 2019).

5 Summary and future directions

In this paper, we presented the mathematical formulation of a new terrestrial biosphere model, QUINCY v1.0, that includes a number of ecophysiological processes (short-term and long-term storage pools, acclimation processes) that have not been represented in earlier TBMs. We evaluated QUINCY v1.0 against a range of observations from worldwide datasets and demonstrated that it is successful in simulating photosynthesis and plant dynamics across large geographical ranges and different ecosystem types with different levels of chronic nutrient input and water availability. We further demonstrated that despite increased complexity and therefore increased numbers of weakly constrained parameters, the model produces predictions of the coupled biogeochemical cycles at site level within reasonable and well-defined bounds. Whether or not the same parameter and therefore process ranking applies also to perturbation situations, e.g. from the manipulation of C, N, or P availability, still remains to be evaluated. The model evaluation provided in this paper points to shortcomings of the baseline QUINCY model in terms of the responses of carbon partitioning to nutrient availability and the interactions between plants and soil processes.

The next logical step in the model's evaluation is to subject it to a range of manipulative experimental settings that will test the realism of the coupling of the carbon–nitrogen–phosphorus cycle in different climate regimes (Zaehle et al., 2014; Yang et al., 2014; Medlyn et al., 2015). Such model evaluation can help point to the adequate representation of individual processes rather than just overall model performance. The inclusion of isotope tracers will increase the scope of such comparisons, since this will allow to better track carbon and nitrogen flows.

One of the motivations behind the development of the QUINCY v1.0 model was the recognition that there may be more than one adequate model representation for one process (Beven, 2008). Therefore, the model has been constructed with a modular design, allowing to replace some of its components by alternative representations to test their effect on model predictions. Examples of future work with QUINCY v1.0 include the control of plants to shape carbon, nitrogen, and phosphorus partitioning and thereby plant functioning in response to environmental change and plant–soil processes by including better representations of the interactions between root growth/activity and SOM turnover and stabilisation.

These basic benchmarks provide a baseline test to integrate QUINCY v1.0 into the JSBACH land surface model (Mauritsen et al., 2019) to allow for a spatially explicit simulation and integration of a range of important processes such as fire disturbance, land-use, and permafrost dynamics. This integration will also allow more comprehensive and rigorous benchmarking against a wider variety of data products. In the midterm, this will allow coupling to the Max Planck Institute's Earth System Model (MPI-ESM) framework ICON to address feedbacks between land biogeochemistry and climate on the basis of an improved representation of biological processes affecting land biogeochemical and biogeophysical processes.

Appendix A: FLUXNET sites

Wohlfahrt et al. (2008b)Beringer et al. (2011)Cleugh et al. (2007)Carrara et al. (2004)Araújo et al. (2002)Dunn et al. (2007)Bergeron et al. (2007)Mkhabela et al. (2009)Ammann et al. (2007)Rebmann et al. (2004)Kutsch et al. (2008)Scherer-Lorenzen et al. (2007)Grünwald and Bernhofer (2007)Lagergren et al. (2008)Sanz et al. (2004)Vargas et al. (2013)Suni et al. (2003)Granier et al. (2000)Berbigier et al. (2001)Keenan et al. (2010)Grünzweig et al. (2003)Tirone et al. (2003)Wohlfahrt et al. (2008a)Tedeschi et al. (2006)Chiesi et al. (2005)Dolman et al. (2002)Lindroth et al. (2008)Lagergren et al. (2008)Goldstein et al. (2000)Urbanski et al. (2007)Hollinger et al. (1999)Schmid et al. (2000)Gu et al. (2006)Scott et al. (2009)Desai et al. (2005)Ma et al. (2007)Ma et al. (2007)Cook et al. (2004)Roupsard et al. (2006)Archibald et al. (2010)

Table A1 FLUXNET sites.