Using Free Air CO₂ Enrichment data to constrain land surface model projections of the terrestrial carbon cycle

Predicting the responses of terrestrial ecosystem carbon to future global change strongly relies on our ability to model accurately the underlying processes at a global scale. However, terrestrial biosphere models representing the carbon and nitrogen cycles and their interactions remain subject to large uncertainties, partly because of unknown or poorly constrained parameters. Parameter estimation is a powerful tool that can be used to optimise these parameters by confronting the model with observations. In this paper, we identify sensitive model parameters from a recent version of the ORgainzing Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) land surface model that includes the nitrogen cycle. These sensitive parameters include ones involved in parameterisations controlling the impact of the nitrogen cycle on the carbon cycle and, in particular, the limitation of photosynthesis due to leaf nitrogen availability. We optimise these ORCHIDEE parameters against carbon flux data collected on sites from the FLUXNET network. However, optimising against present-day observations does not automatically give us confidence in future projections of the model, given that environmental conditions are likely to shift compared to the present day. Manipulation experiments give us a unique look into how the ecosystem may respond to future environmental changes. One such type of manipulation experiment, the Free Air CO₂ Enrichment (FACE) experiment, provides a unique opportunity to assess vegetation response to increasing CO₂ by providing data under ambient and elevated CO₂ conditions. Therefore, to better capture the ecosystem response to increased CO₂, we add the data from two FACE sites to our optimisations, in addition to the FLUXNET data. We use data from both CO₂ conditions of FACE, which allows us to gain extra confidence in the model simulations using this set of parameters. We find that we are able to improve the magnitude of modelled productivity. Although we are unable to correct the interannual variability fully, we start to simulate possible progressive nitrogen limitation at one of the sites. Using an idealised simulation experiment based on increasing atmospheric CO₂ by 1 % yr⁻¹ over 100 years, we find that optimising against only FLUXNET data tends to imply a large fertilisation effect, whereas optimising against FLUXNET and FACE data (with information about nutrient limitation and acclimation of plants) decreases it significantly.

1 Introduction

Since the start of the industrial era, the atmospheric CO₂ concentration has risen from around 278 ppm in 1850 to 417.2 ppm in 2022 (Friedlingstein et al., 2022). Increases in atmospheric CO₂ lead to increases in leaf-scale photosynthesis and intrinsic water-use efficiency, which in turn have the potential to increase plant growth, vegetation biomass and soil organic matter (Walker et al., 2021). Known as CO₂ sequestration, this process transfers carbon (C) from the atmosphere into terrestrial ecosystems. Indeed, terrestrial ecosystems currently remove about 30 % of the CO₂ emitted by human activities each year (Friedlingstein et al., 2020). However, predicting how this carbon sink will evolve under increasing atmospheric CO₂ remains a challenge, especially due to the large uncertainties in the magnitude of carbon-climate feedbacks. Furthermore, the terrestrial ecosystem's ability to store carbon will be influenced by other processes, for example, nutrient limitations (Zaehle and Dalmonech, 2011) – most notably nitrogen (N), which is a key component controlling the carboxylation activity of the RuBisCO in the photosynthetic tissue of the plant.

The large uncertainties in terrestrial carbon projections are largely related to the uncertainty in land surface models, including parametric uncertainty, which relates to the parameter values used in each parameterisation (Zaehle et al., 2005). The first land surface models were developed to provide a physical boundary to meteorology processes. As these models progressed, terrestrial biogeochemical cycles were implemented, simulating leaf gas exchange through Ball–Berry stomatal conductance and plant productivity based on Farquhar photosynthesis (Bonan, 2015). More recently, land surface models have moved from a big leaf model to multi-canopy schemes (Naudts et al., 2015) and started to include the nitrogen cycle and its constraints on the terrestrial carbon balance (e.g. LPJ: Prentice, 2008, OCN: Zaehle and Friend, 2010, ORgainzing Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE)-CN: Vuichard et al., 2019, CLM: Fisher et al., 2019). However, with each new process and complexity added to the model, we add more internal model parameters, which in turn can add more uncertainty. Even though these parameters are generally chosen to represent measurable real-world quantities (e.g. leaf area, plant root depth), their default values are often issued from specific experiments studying the processes at different scales to those used in land surface models. Therefore, it is important to confront simulated model outputs against independent data.

There are a lot of data with which we can evaluate model simulations from vast in situ networks (e.g. FLUXNET, Pastorello et al., 2020) to state-of-the-art satellite retrievals (e.g. Sentinel missions, Malenovskỳ et al., 2012). It is important to evaluate land surface models against these types of data since they help increase confidence in the model simulations. Furthermore, these data can also be used to optimise models through parameter estimation. Parameter estimation methods can be used to perform parameter optimisation where uncertain parameters are tuned to minimise the difference between simulated model output and observed quantities. FLUXNET eddy covariance data have already been used to optimise model parameters in most land surface models, e.g. ORCHIDEE (Kuppel et al., 2012), BETHY (Knorr and Kattge, 2005), JULES (Raoult et al., 2016), Noah (Chaney et al., 2016) and CLM (Post et al., 2017). However, evaluating and optimising against historical trends and present-day observations does not necessarily give us confidence in future projections of the model, given that future environmental conditions are likely to shift compared to the present day (Wieder et al., 2019).

Fortunately, manipulation experiments give us a unique look at how the ecosystem may respond to future environmental change (Van Sundert et al., 2023). One such type of experiment, the Free Air CO₂ Enrichment experiment (FACE; Norby et al., 2010; Walker et al., 2018a), provides a unique opportunity to assess the vegetation response to increasing CO₂. FACE experiments are conducted across several vegetation types and typically consist of two types of plots: one where CO₂ is fumigated to high concentrations and one left as a control.

In particular, 2-decade long FACE experiments in temperate forests of the south-eastern US (Duke and Oak Ridge National Laboratory – ORNL) have been predominately studied to test the representations of carbon–nitrogen cycle processes in land surface models. A full intercomparison of 11 land surface models (Medlyn et al., 2015) demonstrated how these data could be used to evaluate models looking at the effect of ambient and elevated CO₂ on water (De Kauwe et al., 2013), carbon (De Kauwe et al., 2014) and nitrogen (Walker et al., 2014; Zaehle et al., 2014). These two sites were further used in Wieder et al. (2019), where they showed how these experimental manipulations could be incorporated into the model benchmarking tools to help increase confidence in terrestrial carbon cycle projections. FACE experiments can also be used to identify processes that are not well caught by land surface models. For instance, Walker et al. (2019) showed that elevated CO₂ changed carbon allocation to the wood, and none of the models tested was able to reproduce this observation. Combined with warming experiments within a factorial design, the FACE experiments can also be very useful for evaluating how much the models are able to reproduce the single effect of elevated CO₂ vs. the effect of elevated CO₂ when other drivers are changing (De Kauwe et al., 2017). More recently, Sulman et al. (2019) used these two sites to test the effect of adding symbiotic nutrient acquisition strategies to land surface models, and Caldararu et al. (2020) assessed a whole-plant growth optimality approach in improving the representation of leaf nitrogen content compared to existing empirical approaches. The two sites are also the sites we focus on in this study.

We use these sites to check whether the parameterisations and parameters used in a land surface model are able to capture the ecosystem response to increased CO₂. Furthermore, by optimising a land surface model to both ambient and elevated conditions simultaneously, we gain extra confidence in the model projections using this single set of parameters. Although ideally we would want to calibrate under ambient conditions and test the model under elevated conditions, known model structural errors do not guarantee that the model is able to predict changes under different conditions. As such, we provide an alternative approach to model calibration, maximising the available information content of the optimisations. Our study is the first, to our knowledge, to do this with a global land surface model.

Using the ORCHIDEE land surface model as an example, in this paper we show the potential of using manipulation sites to not only optimise unknown model parameters but also increase confidence in the optimised model projections by reducing parameter uncertainty. Furthermore, by optimising parameters linked primarily to the nitrogen cycle as well as considering nitrogen-limited FACE sites, we get an insight into the nitrogen-limiting effect on the fertilising effect of CO₂. This study also looks at how FACE data can complement FLUXNET data in a general optimisation procedure. As such, we aim to answer the following questions.

• Using parameter estimation, can we improve the representation of the simulated productivity of the new nitrogen version of ORCHIDEE over the FLUXNET and FACE sites (under both ambient and elevated conditions)?

• What is the benefit of adding FACE data on top of FLUXNET data when optimising a land surface model?

• How does the future evolution of terrestrial productivity change when simulated using different sets of optimised parameter values?

• Can these experiments help us to describe better the future fertilising effect of CO₂ under possible nitrogen limitation?

2 Methods

2.1 Model

2.1.1 The ORCHIDEE land surface model

The ORCHIDEE model is a global terrestrial ecosystem model developed at the IPSL (Institut Pierre Simon Laplace, France). It simulated the energy (Ducoudré et al., 1993), water (de Rosnay and Polcher, 1998), carbon (Krinner et al., 2005) and nitrogen (Zaehle and Friend, 2010; Vuichard et al., 2019) exchanges between the land surface and the atmosphere. This model can be run at various spatial resolutions, ranging from site to global simulations and over different timescales from 1 d to thousands of years. ORCHIDEE can be run as a stand-alone model driven by meteorological forcing or as part of the IPSL Earth system model (Boucher et al., 2020; Lurton et al., 2020).

In ORCHIDEE, the different types of vegetation are discretised in PFTs (plant functional types) defined by plant metabolism, phenology, type of leaves and local climate. There are a total of 15 PFTs in ORCHIDEE; 8 for the forests, 4 for the grasslands, 2 for the crops and 1 for bare soil. The model describes the different stocks of biomass in the whole soil–plant continuum. There are nine stocks of biomass in the plant: the leaf, the above- and below-ground sapwood, the above- and below-ground heartwood, the fruits, the roots and the long-term and short-term (available to use) reserves. For litter, there are six carbon stocks: metabolic, structural and woody above- and below-ground. Finally, there are four stocks for the soil organic matter: surface, active, slow and passive.

The litter pools are limited by the fall and death of tissues. The pools of organic matter in the soils are alimented by the decomposition of the organic matter in the different pools of the litter. The decomposition of the organic matter is characterised by a fixed residence time for each litter and/or soil pool modulated by environmental conditions.

The $\mathrm{carbon}/\mathrm{nitrogen}$ ratio of leaf biomass is variable, controlled by a supply–demand scheme, while the $\mathrm{C}/\mathrm{N}$ ratio of the other plant pools is a fixed proportion of the leaf $\mathrm{C}/\mathrm{N}$ ratio. A specific $\mathrm{C}/\mathrm{N}$ ratio is set for each soil pool, which varies as a function of the mineral nitrogen in soils. There are also additional mineral nitrogen pools in soils for ammonium, nitrate, nitrous oxides, nitrogen oxides and dinitrogen. The inputs of nitrogen in the soil–plant system are considered to be deposition, fertiliser and manure inputs and biological fixation. Nitrogen losses are associated with leaching, lixiviation and emissions of ammonia, nitrous oxides, nitrogen oxides and dinitrogen.

The nitrogen component for ORCHIDEE was first developed and evaluated inside OCN, a version of the ORCHIDEE model (Zaehle and Friend, 2010). However, at the time, it was not embedded in the operational ORCHIDEE version used in coupled experiments. This component has been recently updated and is now included in default ORCHIDEE runs (Vuichard et al., 2019). This has notably permitted studies of the interactions between the carbon and nitrogen cycles and their effect on gross primary production (GPP). The version of ORCHIDEE we use in our study (ORCHIDEEv3, r6863) is more recent than the one used by Vuichard et al. (2019, r4999). ORCHIDEEv3 (r6863) includes the latest developments of the main ORCHIDEE model (mainly small bug fixes). Furthermore, it includes updates to a few specific N-related processes, notably growth and maintenance respiration. Although this version has been used in the multi-model ensemble for the Global Carbon Budget 2020 (Friedlingstein et al., 2022), it has not yet been optimised against independent data. As such, the initial fit of the model to the FLUXNET data is different than that shown in Vuichard et al. (2019).

2.1.2 Model parameters

An initial list of parameters was compiled based on parameters used in past ORCHIDEE optimisations. This was extended to include parameters of the new nitrogen module selected using expert knowledge of the module developers. Using a Morris sensitivity analysis (Morris, 1991), we remove all parameters to which the different modelled outputs tested (i.e. net primary product NPP and leaf-area index LAI) showed no sensitivity. All the remaining parameters are optimised in this study (Table 1). These parameters represent key parameters of the model-controlling photosynthesis, carbon and nitrogen allocation as well as respiration and global nitrogen cycle behaviour (full descriptions can be found in Appendix A). In addition, the K_Soil parameter is used to control the initial carbon and nitrogen stocks. This parameter makes up for the fact that we cannot reconstruct each site's land-use history and its impacts on the present-day soil carbon stocks. Instead, we add the K_Soil parameter to the optimisation, a multiplication factor applied to some soil carbon and nitrogen pools (slow, passive and labile) to change their initial values. A similar parameter has been used in many previous ORCHIDEE optimisation studies to control the initial carbon stocks of the model (e.g. Santaren et al., 2007; Kuppel et al., 2012; Bastrikov et al., 2018).

For each PFT, the Morris score for each parameter is normalised by the most sensitive parameter. The normalised Morris sensitivity scores are shown in Table 1 and help us understand which are the most sensitive parameters. We see that for sites with a strong seasonal cycle, i.e. TeBS sites, the specific leaf area (SLA) phenology parameters are most sensitive. For the evergreen sites, two of the nitrogen parameters NUE_opt and K_LAtoSA,max gain importance, ranking as highly as SLA.

Table 1 List of parameters used for the optimisation with descriptions, default (prior) model values, ranges of variation and normalised Morris scores denoting the relevant importance of each parameter (labelled “rk” for rank) – darker squares mean “more sensitive”.

2.2 Parameter estimation framework

We perform optimisations by relying on a Bayesian framework to include prior knowledge about the parameters (x_b). Assuming that the errors associated with data observation, model output and parameters follow Gaussian distributions (Santaren et al., 2014), we seek to obtain a posterior optimal parameter set x_opt which corresponds to the minimum of the cost function J(x):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}J\left(\mathit{x}\right)=& \left(M\right(\mathit{x})-\mathit{y}{)}^T{\mathbf{R}}^{-\mathrm{1}}(M\left(\mathit{x}\right)-\mathit{y})\\ +& (\mathit{x}-{\mathit{x}}_{\text{b}}{)}^T{\mathbf{B}}^{-\mathrm{1}}(\mathit{x}-{\mathit{x}}_{\text{b}}).\end{array}\end{array}$$

For a given parameter set x, J(x) measures the mismatch between observations y and the corresponding model outputs M(x) as well as the mismatch between the prior, or background, parameter set x_b and x. Each of these terms is weighted by its error covariance matrices, R and B, for the observations and parameters, respectively (Tarantola, 2005). In this study, we set both matrices to be diagonal. For B, we define the prior distribution of each parameter to be 40 % of the prior range. For R, we define the observation error (variance) as the mean-squared difference between the observations and the prior model simulation so that this variance reflects not only the measurement errors, but also the model errors. Furthermore, since we do not consider error covariances, R is diagonal, and therefore we can decompose the first term of Eq. (1) into different terms for each assimilated data stream:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}J\left(\mathit{x}\right)=& K_{\text{Flx}}\left(M_{\text{Flx}}\right(\mathit{x}\left)\mathrm{-}{\mathit{y}}_{\text{Flx}}{)}^T\right({\mathit{\sigma }}_{\text{Flx}}^{-\mathrm{1}}\left)\right(M_{\text{Flx}}\left(\mathit{x}\right)\\ -& {\mathit{y}}_{\text{Flx}})+K_{\text{FACE}}(M_{\text{FACE}}\left(\mathit{x}\right)\\ \mathrm{-}& {\mathit{y}}_{\text{FACE}}{)}^T\left({\mathit{\sigma }}_{\text{FACE}}^{-\mathrm{1}}\right)\left(M_{\text{FACE}}\right(\mathit{x})\\ -& {\mathit{y}}_{\text{FACE}})+(\mathit{x}-{\mathit{x}}_{\text{b}}{)}^T{\mathbf{B}}^{-\mathrm{1}}(\mathit{x}-{\mathit{x}}_{\text{b}}).\end{array}\end{array}$$

The Flx and FACE subscripts are used to denote the FLUXNET and FACE parts of the equation; K_i denotes the weighting using for each data stream, σ_i denotes the observational error, and M_i and y_i denote modelled and observed data streams.

There exist many different approaches we can use to find the set of parameters which minimise J(x). These range from simple manual tuning, which are very computationally demanding and inefficient, to more complex algorithms either based on deterministic gradient descent methods or stochastic random search methods. Using “ORCHIDAS”, the ORCHIDEE data assimilation tool developed at the Laboratoire des Sciences du Climat et de l'Environnement (Bastrikov et al., 2018), we performed a couple of preliminary experiments to determine which algorithm to use. We tested a gradient descent method based on the L-BFGS-B algorithm (limited memory Broyden–Fletcher–Goldfarb–Shanno algorithm with bound constraints BFGS; Byrd et al., 1995) and a random search method based on the genetic algorithm (GA; Goldberg and Holland, 1988; Haupt and Haupt, 2004). We found that the GA method outperformed the gradient method in reducing the cost function. These initial results are coherent with the Bastrikov et al. (2018) study, which optimised the GPP and latent heat fluxes of a former version of ORCHIDEE against a number of FLUXNET site measurements and also found that the GA algorithm outperformed the other methods, notably by allowing a full exploration of all the parameter space.

The genetic algorithm consists in applying the laws of evolution to our set of parameters by considering the set of parameters as a chromosome, with each parameter as a gene. At each iteration, the algorithm fills h chromosomes with parameter values. The first pool of chromosomes is created by randomly perturbing the value of the parameter. For the following iterations, the chromosomes are created from the previous iterations' chromosomes. Two processes come into play: (a) a crossover process, where we have an exchange of genes between two chromosomes; and (b) a mutation process, where random genes are perturbed. To ensure that the best chromosomes get the most descendants, each chromosome of each iteration is ranked as a function of the cost associated with the parameter's value in the chromosome.

2.3 In situ data

In this study, we consider two sites from the FACE network in nitrogen-limited temperate forest ecosystems: Oak Ridge (ORNL; Norby et al., 2010), a site dominated by broadleaf deciduous forests (TeBS, for temperate broadleaf summergreen forests); and Duke (DUKE; McCarthy et al., 2010), a site dominated by needleleaf evergreen forests (TeNE, for temperate needleleaf evergreen forests). The data for these sites come from the FACE Model Data Synthesis project (Walker et al., 2018a, b, https://facedata.ornl.gov/facemds/, last access: 15 February 2024). For each site, we use the data from two experimental plots (with their associated error bars), one with unaffected atmospheric CO₂, i.e. ambient (AMB), and one with elevated atmospheric CO₂ (ELE). Although the DUKE experiment also has ammonium nitrate treatments at half of its plots from 2005 onwards (Feng et al., 2010), we only consider the data from the plots without nitrogen fertilisation.

The version of ORCHIDEE we use in this study has yet to be optimised against FLUXNET data using a Bayesian framework, as was done with previous nitrogen-free versions of the model (e.g. Kuppel et al., 2012; Peylin et al., 2016). Therefore, we also consider TeBS and TeNE sites from the FLUXNET2015 dataset (Pastorello et al., 2020). This dataset provides gap-filled half-hourly meteorological data measured at each site (air temperature, humidity, pressure, wind speed, rainfall and snowfall rates, short-wave and long-wave incoming radiation; see Vuichard and Papale, 2015). It also provides net carbon flux measurements, such as net ecosystem exchange (NEE) further split into GPP and total ecosystem respiration (TER) following a classical night-time vs. day-time flux partition (Lasslop et al., 2010). For each of the two vegetation types, sites with over 60 % vegetation coverage are kept. We exclude sites with overly large discrepancies with the prior model output, such as with no apparent seasonal cycle, large data gaps or only 1 year of data. The list of in situ sites used can be seen in Table 2 partitioned by vegetation type.

Table 2 List of the in situ FACE and FLUXNET sites used in the study. The FLUXNET sites are labelled by a country code (first two letters) and site name (last three letters). The FACE sites are both found in the US. The period corresponds to the available years of data for each of the sites.

2.4 Performed experiments

Before performing the optimisations, for each of the sites in this study, a two-step spin-up is performed. The first step helps to put the prognostic variables, including vegetation state, soil carbon pools and soil moisture content, at equilibrium. The available meteorological forcing is cycled over several millennia (with pre-industrial CO₂ concentrations) to ensure that the long-term net carbon flux is close to zero. After reaching equilibrium, a second simulation is performed (transient) from the year 1860 to 1 year before the first forcing year while increasing the CO₂ concentration at each simulation year following global historical observations.

Before performing the optimisations, we also conduct a sensitivity analysis of the parameters (as described in Sect. 2.1.2 and shown in Table 1). A sensitivity analysis tests how differently the model outputs change with respect to different parameters. This is done to ensure that only parameters showing some sensitivity to the model outputs are used in the optimisation, thereby minimising the risk of using parameters that are weakly constrained by the fluxes. This is an important step since we want to avoid constraining parameters that will have a small impact on the optimisation but that have the potential to significantly degrade the model–data fit against processes not included in the calibration.

Once spun up and with the list-sensitive parameters, we perform two main sets of optimisations, always starting from this spin-up. The first is over the FLUXNET sites only, while the second also includes data from the FACE sites. Due to the CO₂ fumigation over the FACE sites, NEE is not measured at these sites, and therefore GPP and TER estimates cannot be derived. Instead, for the FACE sites, we have annual NPP and daily LAI data. Throughout this study, we perform multi-site (MS) optimisations, i.e. optimisations executed over multiple sites of the same PFT simultaneously in order to find one common set of optimised parameters. Each optimisation is run for 20 iterations, which we found to be sufficient for the system to converge. For each iteration, 32 chromosomes are used, i.e. 32 different combinations of parameter values. We leave the last year of each FLUXNET site out of the optimisation to have independent data for the validation step of the analysis.

The first set of optimisations tests two different combinations of gross and net carbon fluxes.

• Flx_GR: two MS optimisations against daily GPP and TER, one for all of the TeBS sites and one for all of the TeNE sites

• Flx_GN: two MS optimisations against daily GPP and NEE, one for all of the TeBS sites and one for all of the TeNE sites

In each case, two fluxes are used in optimisations. Note that GPP and TER are derived from NEE with NEE = TER − GPP. This means that these data are model-derived estimates, which could introduce additional uncertainty into the results. However, by separating the fluxes, we get a better understanding of the underlying mechanisms constraining two different ecosystem functions and are able better to diagnose the overestimations or underestimations of the assimilated processes, as initially discussed in Santaren et al. (2007). We are especially interested in the GPP constraint since this will give us an insight into plant productivity and will allow us to assess the CO₂ fertilising effect under nitrogen limitation. GPP is also directly used in the calculation of water use efficiency (WUE), here defined as the ratio between GPP and transpiration, one of the diagnostics we consider at the end of the study.

We further acknowledge that the data streams are not independent of one other. This poses a challenge when working in a Bayesian framework, especially when defining the R matrix in Eq. (1). Although there are methods for including the correlation between different data streams in R, these are relatively new and require a lot of extra analysis beyond the scope of this study. Instead, we rely on the standard method of inflating variances (Chevallier, 2007).

The optimal parameters found by optimising against the FLUXNET sites improve the fit to contemporary data. However, it is unclear whether the predicative skill of the model is improved. Therefore, after assessing the FLUXNET results, the next step is to incorporate the FACE sites. Using a simultaneous approach, the FACE and FLUXNET sites are optimised together in this second set of experiments. This approach ensures that the information is not lost between steps, as could be the case in the step-wise approach when the optimisations are done one after the other. The optimisations are set up to give a higher weight to the single FACE site in each case, so that K_Flx=1 and K_FACE=n in Eq. (2), where n is the number of FLUXNET sites for the given PFT. Based on our results (see Sect. 3.1) and our motivation to better capture the productivity of the different ecosystems, we choose to focus on the former FLUXNET optimisation, i.e. the one against GPP and NEE. Each of the following FACE site experiments is performed simultaneously with a Flx_GN optimisation over the relevant PFT.

• Flx_GN-AMB: two optimisations against annual NPP and daily LAI, one each for the DUKE and ORNL sites at ambient CO₂ concentrations, perform simultaneously with a GPP–NEE multi-site FLUXNET optimisation.

• Flx_GN-ELE: two optimisations against annual NPP and daily LAI, one each for the DUKE and ORNL sites at elevated CO₂ concentrations, perform simultaneously with a GPP–NEE multi-site FLUXNET optimisation.

• Flx_GN-BOTH: two optimisations against annual NPP and daily LAI, one each for the DUKE and ORNL sites with both ambient and elevated CO₂ concentrations simultaneous, perform simultaneously with a GPP–NEE multi-site FLUXNET optimisation.

For the final part of this study, we consider the sensitivity of the simulated GPP, NPP and WUE to CO₂ increase whilst keeping the other drivers constant. Each of the FLUXNET sites is tested by running idealised 100-year long simulations starting from present-day atmospheric CO₂ and 380 ppm and increasing CO₂ by 1 % yr⁻¹, leading to a near tripling of CO₂ by the end of the simulation. This is done for both the prior and optimised models, using default model parameters and the Flx_GN-BOTH model parameters, respectively.

3 Results and discussion

3.1 FLUXNET optimisations

The nitrogen cycle version of the ORCHIDEE model used in this study has not yet been optimised against FLUXNET data, although it has been extensively tuned manually. Therefore, the first step is to see whether the fluxes over the TeBS and TeNE sites are well represented in the model and whether they can be optimised using observed data and the sensitive parameters identified in the Morris experiment (Table 1). For the FLUXNET optimisation, two tests are conducted. The best-performing optimisation will serve as the starting point for the optimisations including data from the FACE sites. In this section, we present the results from both MS FLUXNET optimisations: Flx_GR and Flx_GN. Figure 1 shows the mean seasonal cycle across all the FLUXNET sites for each PFT considered. We show the modelled GPP, TER and NEE fluxes against the observed time series.

Figure 1 The main panel in each column shows the PFT-averaged mean seasonal cycles of daily observed and simulated GPP, TER and NEE fluxes using different parameter values. The side panels show the model–data RMSD for the daily time series at each site, with the black horizontal bars showing the mean value across the sites.

For the deciduous sites (TeBS; Fig. 1, left-hand column), we see that both GPP and TER are overestimated by the prior model, and the NEE sink is underestimated. This overestimation is the most severe for TER, where the prior model simulates a maximum of approximately 9 gC m⁻² d⁻¹ when the maximum TER observed is half that. In contrast, the overestimation for GPP is very slight when found at and after the peak. Both optimisations improve the model–data fit against GPP by correcting the overestimation found after the peak. Flx_GN performs best, with the average seasonal peak now the same as the observations. Flx_GR on the other hand reduces the peak below that observed. Similarly, both optimisations now start production later in the year, degrading the fit to the observations in the early months. In ORCHIDEE, deciduous sites lose all their leaves in winter, and therefore no photosynthesis occurs before the leaves start growing back in spring. In contrast, the observations never go to zero, implying that there is undergrowth or evergreen vegetation present that we are not accounting for in the model set-up. When looking at the root-mean-square deviations (RMSDs) of the individual sites, we also see that Flx_GN reduces the spread relative to the prior. For TER, both optimisations improve the model–data fit over the whole period, with Flx_GR performing slightly better. This is not surprising since this optimisation directly considers the TER component of NEE. The optimisations mainly change the magnitude of the peak and do not correct for its late timing. When looking at the RMSD, we can see that a group of sites is driving the overestimation of TER, with values close to 4.5 – this is corrected for in both optimisations. For NEE, the Flx_GN optimisation performs better than Flx_GR, especially in fitting the autumn month. However, the overestimation of summer TER with this parameter set means we do not attain the minimum of the NEE trough. We note that, for NEE, the posterior spreads of the RMSDs over all the sites are the same for both optimisations.

For the evergreen sites (TeNE; Fig. 1, right-hand column), the optimisations improve the GPP underestimation in the prior model by increasing the peak by approximately 2.5 gC m⁻² d⁻¹. However, this falls short of correcting the full overestimation, which is closer to 4 gC m⁻² d⁻¹. The Flx_GN optimisation performs best, with the timing of the average peak closest to the observed value. However, both optimisations degrade the average fit to TER, increasing the overestimation found when using the prior set of parameters. Both optimisations move the summer peak up by between 0.5 and 1 gC m⁻² d⁻¹, with Flx_GN increasing the most. When considering the fit of the individual sites (right-hard part of each panel), we note that two anomalous sites are driving this behaviour. These sites (IT-Lav and US-Wi4) have respiration rates much lower than the other sites. Since we cannot match the respiration rates, we cannot capture the full NEE dip in summer. However, the improved GPP with both optimisations does mean that the NEE is slightly improved compared to the prior mode.

In these optimisations, we have included a parameter, K_Soil, which acts as a multiplicative factor of the initial soil and nitrogen pools (slow, passive and labile). We have used such a parameter in past ORCHIDEE optimisations at the FLUXNET sites (e.g. Kuppel et al., 2012; Peylin et al., 2016; Bastrikov et al., 2018; Bacour et al., 2023) and found that it played a large role in improving the model–data fit against respiration. Therefore, it seems counter-intuitive that we do not improve the fit to respiration as much as expected when including K_Soil in the ORCHIDEE v3 optimisations, especially for the TeNE sites. The past ORCHIDEE experiments all used a previous version of ORCHIDEE without the nitrogen cycle, so this factor solely acted on the carbon pools and heterotrophic respiration. Here K_Soil multiplies both carbon and nitrogen pools to maintain the $\mathrm{carbon}/\mathrm{nitrogen}$ ratio. However, by acting on the nitrogen pools, we directly impact the mineralisation rate and thus indirectly the plant N uptake, leaf N content, V_max and, therefore, GPP. While K_Soil used only to impact soil respiration, it now impacts both respiration and GPP, and so the optimisation needs to find a compromise to fit both. To adjust the respiration, K_Soil is decreased, reducing the carbon and nitrogen pools in the soils, but, at the same time, GPP is significantly reduced, deteriorating its fit to observations.

Overall, the ORCHIDEE model reasonably represents the TeBS and TeNE carbon fluxes, although respiration at the TeNE sites is high, even after optimisation. The Flx_GN optimisation results in the best-simulated production for both types of vegetation.

3.2 Incorporating data from the FACE sites

Given the results from the previous section and our motivation to improve the model performance regarding ecosystem productivity, we will further include the FACE data in the Flx_GN optimisation.

Figure 2 Time series of annual NPP under ambient (NPP_AMB; top) and elevated (NPP_ELE; middle) CO₂ conditions for each of the ORNL (a) and Duke (b) FACE sites. The ratios between NPP under elevated conditions and NPP under ambient conditions are shown in the bottom row, with a dotted-grey line indicating where the ratio is 1. The observations are shown in black with error bars. The coloured lines represent different ORCHIDEE model simulations under different parameter sets – “prior” refers to the standard ORCHIDEE run using uncalibrated parameters, and the different “Flx” runs denote the data streams used in the optimisations (defined in Sect. 2.4). The bar chart inset in each time series panel shows the RMSD between the observations and each model run. The multi-annual NPP mean ± its standard deviation is shown next to each panel.

3.2.1 Improving simulated NPP values

At ORNL, we can see in Fig. 2 that the uncalibrated version of ORCHIDEE (prior) overestimates the yearly NPP under both ambient and elevated conditions. This is consistent with prior GPP overestimation observed at the TeBS FLUXNET sites (Fig. 1). When using parameters from the FLUXNET-only optimisation (Flx_GN), we partly reduce this overestimation. For both CO₂ conditions, including FACE data as an additional constraint in the optimisation (Flx_GN-AMB, Flx_GN-ELE and Flx_GN-BOTH) improves the estimation of NPP compared to solely relying on the data from the FLUXNET sites. Under all the atmospheric CO₂ conditions, Flx_GN-ELE reduces the RMSD the most, followed by Flx_GN-BOTH and then by Flx_GN-AMB. We would expect the latter to perform best at fitting NPP_AMB since it uses the observations in the optimisation, unlike Flx_GN-ELE. However, as we will see later in Fig. 3, this is because the fit to LAI, the other part of the cost function, is improved most with the Flx_GN-AMB optimisation.

Figure 3 Model fit against seasonal LAI under different parameter sets at ORNl and DUKE. LAI under ambient CO₂ conditions (LAI_AMB) shown in the top row and elevated CO₂ conditions (LAI_ELE) in the bottom row. The bar chart inset in each panel shows the RMSD between the observations and the model.

When considering the NPP ratio (elevated over ambient values) at ORNL, we see that the observations show a decreasing trend, suggesting a possible progressive nitrogen limitation. The prior model is unable to capture this trend with a fixed ratio of around 1.3. Similarly, the FLUXNET-only and optimisation under ambient conditions do not mimic this decreasing trend. When using the FACE data from the elevated experiments in optimisation (Flx_GN-ELE and Flx_GN-BOTH), the resulting ratio trend is negative. The observations of NPP_AMB and NPP_ELE clearly show a slight decrease in NPP for the last years of the observations. Based on additional experiments, it has been shown that this limitation is due to nitrogen limitation (Norby et al., 2010). Again, the prior does not capture the NPP decrease. For optimisations, Flx_GN-ELE and Flx_GN-BOTH may show a small signal, although it is not so obvious. This leads us to think that ORCHIDEE does not reproduce the dynamics of nitrogen in the soil well, notably the reduction in nitrogen availability for the plant, at least for this site. It is also possible that the model starts with too large a nitrogen content in the soil. Since parameter optimisation is insufficient for capturing this trend, we would instead need structural changes to the ORCHIDEE land surface model or possibly set a better initial nitrogen available content.

Figure 4 Posterior parameter values of the FLUXNET-only optimisation (blue, Flx_GN) and the optimisation incorporating FACE data under both atmospheric regimes (orange, Flx_GN-BOTH). The horizontal line represents the prior value of each parameter and each box spans the range of variation allowed for each parameter during the optimisation. The parameters highlighted in bold correspond to the parameters identified as most sensitive in the preliminary Morris experiment (see Table 1).

For DUKE, the prior model underestimates the NPP values under both CO₂ conditions, which again is consistent with the GPP fit seen at the other TeNE sites (Fig. 1). The optimisation against FLUXNET data corrects this underestimation slightly. However, it is only by including the FACE data in the optimisation that we get the right magnitude and start getting the right inter-annual pattern. In each case, the optimisation performing best is the one that includes the relevant observations, i.e. Flx_GN-AMB under ambient conditions and Flx_GN-ELE under elevated conditions. The Flx_GN-BOTH optimisation is able to fit both the ambient and elevated data, reducing the RMSD to a similar extent to the best optimisation. When considering the elevated-to-ambient NPP ratio for DUKE, we see that the prior model and the FLUXNET-only optimisation (Flx_GN) perform worst, with some values below 1. Values below 1 suggest that the forest is less productive under increased atmospheric CO₂. Looking more closely at the ORCHIDEE simulations for these years and parameter settings, we found that the mean-annual GPP under elevated conditions is higher than GPP under ambient conditions. Therefore, changes in autotrophic respiration are responsible for the negative NPP ratio at the site. Maintenance respiration is a function of leaf nitrogen content; the more nitrogen, the more respiration. It is, therefore, possible that the maintenance respiration sensitivity to leaf nitrogen is too high in the model for TeNE under these parameter settings, at least at the DUKE site. The optimisations using FACE data do much better at simulating positive ratios, with Flx_GN-BOTH performing best. However, none of the optimisations constantly achieves the high magnitude of around 1.3 in the observations.

Figure 5 Fractional change in the model–data RMSD for the daily NEE and GPP at each site grouped by PFT obtained by the different optimisations. Fractional change is calculated by dividing the posterior RMSD by the prior RMSD (i.e. RMSD$\left[y,M\left(x_{\text{opt}}\right)\right]/\text{RMSD}\left[y,M\left(x_{\text{b}}\right)\right]$ using the notation in Sect. 2.2). Values below 1 represent an improvement in the model–data fit, whereas values above 1 represent an increase in RMSD compared to the prior run. The crosses show the mean value across the sites for a given optimal parameter set. This figure complements Fig. 1, here showing the optimisations using the FACE data fit the time series compared to the prior model and the optimisation (GN) without FACE.

3.2.2 Improving the fit of the LAI

During the FACE experiments, we also optimised the model against daily LAI values (Fig. 3) since this allows us to capture some of the sub-annual variability driving NPP. For ORNL, under ambient and elevated CO₂, the prior model overestimates the LAI, peaking late and keeping its leaves late into the winter. The optimisation against the FLUXNET sites partly fixes both of these issues, bringing forward the leaf fall and decreasing the seasonal peak to 65 % of the observed peak. When adding the FACE data under ambient conditions to the optimisation (i.e. Flx_GN-AMB) we get the best fit, although the peak is still underestimated. The optimisation against both atmospheric conditions (Flx_GN-BOTH) performs similarly to Flx_GN-AMB. When we only add the elevated data to the optimisations, we do the worst of the FACE optimisations, with the summer peak nearly half of the observed one. This explains why Flx_GN-ELE outperformed Flx_GN-AMB in fitting NPP at ORNL (Fig. 2). By improving NPP to a larger extent than Flx_GN-AMB, Flx_GN-ELE was unable to improve the LAI as much. This highlights one of the common issues with multi-stream parameter estimation – improving the fit against one data stream can be part of one trade-off against another. Sometimes this can even be a degradation compared to the prior. Finally, although leaves are kept into the winter, unlike in the observations, the LAI seasonality is much improved in the optimisations compared to the prior model.

3.2.3 Posterior parameter values

In Fig. 4, we consider how the different parameters have changed by comparing results from the FLUXNET-only optimisation and the optimisation including the FACE data. For the temperate broadleaf summergreen parameters it is hard to distinguish significant patterns. The posterior values sometimes move in the same direction for both optimisations (e.g. K_N, FCN_wood, R_root and z), and at other times the parameters are pulled in different directions (e.g. CN_leaf,max and R_leaf). When considering the most sensitive parameters, both optimisations agree with the direction of change, with the exception of NUEopt. For the photosynthetic parameters A₁ and B₁, we see that these are changed during the FLUXNET-only optimisation but are less important in correcting the fit when the FACE data are included in the optimisation. Overall, this highlights that adding the FACE data can pull the parameters in a different direction than when using only FLUXNET data. This indicates that there are likely compensating effects (or a different repartitioning of the model error). It is possible that, with more datasets, the results would be more robust, but this would need to be confirmed with the use and assimilation of additional data streams.

In contrast, for the temperate needleleaf parameters, nearly half (11 out of 23) of the parameters were unchanged (or only slightly changed) during the FLUXNET-only optimisation. However, when optimised with the additional FACE data, these parameters changed greatly. These parameters include a number of the most sensitive parameters, suggesting that these are especially important in capturing the model response under both atmospheric regimes for needleleaf sites.

Although looking at these parameter values can be very informative, we must remember that there are complex interactions between the parameters and processes that will not be evident from just looking at these values.

3.2.4 Maintaining the fit to the FLUXNET sites

The misfit part of the cost function (i.e. measuring the difference between the model and observations) is made up of two components during the FACE optimisations. The first is the fit of the FACE site to annual NPP (Fig. 2) and daily LAI (Fig. 3) under ambient or evaluated conditions. The second is the fit of the multiple FLUXNET sites to daily GPP and NEE. Although the FACE component was weighted during the optimisation (see Sect. 2.4), it is important that we maintain a good model–data fit for the FLUXNET sites to ensure confidence in the optimised parameter values. In Fig. 5, we consider the fractional change in the model–data RMSD (calculated by dividing the posterior RMSD by the prior RMSD). In other words, this helps us see to what extent the calibration with FACE data changes the estimated fluxes shown in Fig. 1. For both vegetation types, the best improvement against both GPP and NEE is found when using parameters from the experiment solely optimising over the FLUXNET sites (Flx_GN). This is true for both the optimisation years, i.e. the years used in the optimisation, and the validation year, i.e. the year of independent data omitted from the optimisation. Adding further constraints by including FACE data reduces the improvement in the model–data fit. After all, even though they are closely related, the FACE optimisations looked at improving the fit to NPP and LAI, not GPP and NEE. However, the sets of parameters obtained from these latter optimisations will give a better compromise for both the FACE and FLUXNET sites.

For the TeBS sites, on average, we improve the fit to GPP over the optimisation years to the same extent for all the optimisations. Slightly higher reductions in the fractional GPP RMSD change are observed for the validation year. However, for NEE, fractions in RMSD are different depending on the optimisation – with the ambient and both optimisations resulting in the smallest decreases for the calibration years. For the validation, these two optimisations degraded the fit for the validation year (20 %).

In contrast, for the TeNE sites against the calibration years, the optimisations all give similar reductions in the model–data fit RMSD for both NEE and GPP. When confronted with the validation year, the response is more spread out, with the FLUXNET-only optimisation performing best, closely followed by Flx_GN-BOTH.

For the validation, we only used 1 year of data, which helps to explain why the results are more spread out. With 1 year of data, we are more susceptible to specifics of that given year instead of trends over a longer period. Ideally, we would want to validate over a larger period. However, some of the sites in this study only had a few years in their data record, making this impossible.

3.3 Projections using the optimised models

To conclude this study, we test how the optimised model parameters impact the model responses to CO₂ increase. We especially want to consider the additional information gained from including FACE data in the optimisation. As such, for this last experiment, we use parameters of the Flx_GN (FLUXNET-only) and Flx_GN-BOTH (FLUXNET and FACE data) optimisations. The Flx_GN-BOTH optimisation resulted in the best compromise in simulating NPP under ambient and elevated conditions as well as their ratio while also maintaining a satisfying model–data fit to GPP over the FLUXNET sites.

Table 3 Mean values at the start of the 100-year experiment and net change by the end of the experiment for the prior model and the optimised models.

Figure 6 Effect of changes in the atmospheric CO₂ concentration on GPP, NPP, water use efficiency (WUE) and plant N uptake for different model parameter sets. Atmospheric changes are represented by a 1 % CO₂ increase per year. The prior model (grey) and optimised model (Flx_GN-BOTH; orange) were run at each FLUXNET site (see Table 2). The thick lines represent the mean simulated across all the sites, while the shaded areas represent the standard deviation. The mean and standard deviation at the end of the 100-year simulation are shown on the right-hand side of each panel.

To test the model response to CO₂ increase, we run an idealised 100-year long experiment increasing CO₂ by 1 % yr⁻¹ over all the FLUXNET sites. This is similar to the experiment done in Vuichard et al. (2019), where the model was also assessed with GPP, transpiration and WUE. Note that the prior projections shown here are similar to those obtained in Vuichard et al. (2019). We first note that the optimised models predict lower starting values of GPP, transpiration and WUE for TeBS when compared to the prior model (Table 3). This is consistent with Figs. 1 and 2, where for TeBS the prior had the highest productivity. Similarly, the behaviour for the TeNE sites mirrors earlier findings – the prior model underestimated productivity. The optimisations lead to higher starting values (except for WUE), and including data from FACE is found to induce a larger increase. When considering the total plant N uptake for both vegetation types, the prior starts with a high value and the optimised runs with much lower values.

In Fig. 6, we consider the net increase over the 100-year simulation. For both vegetation types, the rate of increase in GPP found by using parameters from the FLUXNET-only optimisation is the highest. For TeNE, the curve is starting to plateau due to a progressive nitrogen limitation. However, in both cases, the rate of increase is much faster compared to the runs with parameters from the Flx_GN-BOTH optimisation. The fertilisation rate decreases compared to the FLUXNET-only optimisation. In Fig. 2, we saw that the Flx_GN-BOTH optimisation was starting to model the declining productivity response at ORNL. We may still be overestimating the fertilisation effect over the TeBS sites, and we could expect an even stronger response to nitrogen limitation – further reducing this GPP increase. Although the GPP increase is similar to the prior model for forest types, the larger spread among the sites for Flx_GN-BOTH is notable – with a negative GPP change at some sites suggesting a stronger nitrogen limitation under these parameters.

For transpiration rates over both TeBS and TeNE, the Flx_GN simulations result in the smallest change by the end of the century. In contrast, the simulations Flx_GN-BOTH result in the largest change. For WUE, the differences observed in GPP and transpiration are mostly cancelled out, and therefore the different parameters do not elicit a great variation in responses. For TeBS, the prior and Flx_GN simulations give a similar increase, with Flx_GN-BOTH suggesting a slightly weaker increase. For TeNE, it is the opposite – Flx_GN is slightly weaker than the other two simulations.

For plant N uptake, changes over the 100-year period are small in magnitude but vary between the different optimisations. The prior model shows a range of responses over the different sites (evidenced by the large spread), overall decreasing the rate of N uptake for TeBS and increasing the rate of uptake for TeNE. The FLUXNET-only optimisations lead to a slight increase in plant N uptake by the end of the century for both vegetation types. In contrast, the FACE optimisation leads to a decreased plant N uptake. This is especially notable for TeNE, where the spread of responses over the sites is greatly reduced. In all the cases, changes to the rate of change occur during the first half of the simulations, plateauing to a constant value for the rest of the runs.

Note that our experiment only looks at increasing CO₂ while keeping the other drivers constant. It is possible that we would see different responses if we were to include changes in meteorological forcing (to mimic climate change) or changes in nitrogen deposition and fertilisation, which would change nitrogen limitation and responses to water stress. Although we performed a sensitivity analysis to select the parameters, it is possible that additional parameters (e.g. ones controlling water stress or controlling the allocation of carbon in the plant) will give a different response. Furthermore, direct structural changes or addition to the model code could also change the results. Nevertheless, the fact that different parameters give such a varied model response to elevated CO₂ for a given model structure highlights the importance of finding a robust set of parameters to have faith in. When performing the FLUXNET optimisations, the optimisation against GPP and NEE gave the best fit to both quantities, resulting in sets of parameters that worked well across each PFT. However, when considering production through NPP and LAI under different atmospheric CO₂ conditions, we found that the parameters were unable to capture the differences observed in data under different CO₂ regimes. Instead, the parameter sets found by optimising the model against FLUXNET data and both ambient and elevated CO₂ conditions were found to be more robust. These parameters resulted in the best fits overall against NPP and LAI under two different CO₂ regimes. As such, we have more confidence in these parameters and their ability to simulate terrestrial production under different atmospheric CO₂ conditions, leading us to have more faith in the model projections performed with them.

4 Conclusions

As our terrestrial models become more complex through the addition of more processes, we need to confront them with observed data to ensure we have confidence in the model's predictions. Manipulation experiments allow us to test the model under different CO₂ regimes and its capabilities to reproduce the ecosystem responses. FACE sites in particular are an important tool in evaluating modelled ecosystem response to climate change. They can be thought of as space-for-time substitution experiments (Rastetter, 1996) but where the change in atmospheric CO₂ is controlled and even manipulated to exceed conditions naturally found around the globe currently. By optimising model parameters against data from both ambient and elevated atmospheric CO₂ conditions, we gain confidence in the model's ability to reproduce fluxes under different atmospheric conditions. It will be interesting to use these parameters further to assess the evaluation of carbon stocks under high concentrations of CO₂ and at a larger scale to evaluate more directly the impact on the global sink.

We find that through the different optimisations of this carbon–nitrogen version of ORCHIDEE, we are able to improve the representation of simulated productivity. All the optimisations are able to improve modelled GPP, and we generally improve the magnitude for NPP for the two levels of atmospheric CO₂. However, we do not achieve as good an improvement against respiration and, therefore, against NEE. Although we are unable to capture fully the inter-annual variability of NPP after optimisation, at ORNL we do start to model a negative trend for the NPP ratio, which is apparent in the observations but is not simulated in the prior model. This suggests that the optimised parameters are able to capture the progressive nitrogen limitation at this site. We do struggle to capture the seasonal cycle of LAI properly at both FACE sites, suggesting incorrect LAI allocation. These results highlight the fact that optimising a land surface model (LSM) with the nitrogen cycle is more difficult and complex than with a carbon-only LSM given the increased model feedbacks. In particular, the dependence of plant productivity on soil nitrogen availability and in contrast the dependence of soil N content on litter input (and hence productivity) induce strong positive feedbacks. This provides a warning to other modelling groups looking to calibrate the carbon and nitrogen cycles in their models. Overall, the current optimisation performs slightly worse compared to the prior optimisations of previous ORCHIDEE versions (e.g. Kuppel et al., 2012). One example of increased model feedbacks is through the K_Soil parameter. Without the nitrogen cycle, by changing the initial carbon stocks, this parameter was able to fix the magnitude of the respiration flux. However, the parameter now also changes the initial nitrogen stock and hence the mineralisation flux in the soil, which impacts GPP. Another approach would have been to have several multiplicative factors, each changing different pools, or indeed one for each C and N. However, this would likely lead to more complications given the strong feedbacks observed. If one pool declines more than others in terms of C and N content or if one pool becomes more depleted in N, it is probable that the model will enter a transient phase with potentially strong compensating fluxes, such as a large net carbon flux, to restore internal consistency. Ideally, we would optimise various model parameters that govern the turnover time and the $\mathrm{C}/\mathrm{N}$ ratio of each pool throughout the entire spin-up period. However, achieving such optimisation is not currently feasible.

Although the optimisation is not as optimal as that achieved with a carbon-only model, this work opens a new avenue to validate LSMs quantitatively with FACE data. We see that not only is there a benefit to adding FACE data on top of FLUXNET data when optimising a land surface model, but that it is, in fact, risky not to. The FLUXNET-only optimisations do not perform well under elevated conditions, which is critical when predicting the terrestrial response to climate change. Furthermore, since we see that the future evolution of terrestrial productivity change is sensitive to the parameter values used in the model, getting these parameters right is critical. This is notable for both vegetation types where the FLUXNET-only optimisations and the optimisations with FACE data give very different trajectories in the idealised 1 % CO₂ experiments, with the FLUXNET-only optimisations likely overestimating the CO₂ fertilisation effect. However, we do need to be cautious in assessing these results since we are only using one FACE site for each PFT, meaning we are likely tuning to the specificities of that site. For example, ORNL shows a progressive nitrogen limitation, but this is not expected over all the sites. Ideally, we would include a lot more FACE sites to capture different conditions. In particular, if we could optimise by grouping sites based on different levels of nitrogen limitation, then if the posterior parameters are found to be similar, the model processes allow for these differences.

In any optimisation, there is always a danger of overfitting to data, limiting the generalisability of the calibrated model. By optimising the model against a number of different constraints (i.e. more than one data stream), we decrease the risk of overfitting and therefore gain some confidence in our parameter and hence in the projections. Such experiments can help us to describe better the future fertilising effect of CO₂ under possible nitrogen limitation. However, we find that, in our study, two sites are not sufficient to draw such conclusions about terrestrial responses to elevated CO₂, which could vary over different ecosystems. Although we have shown this approach of joint optimisations to be promising, more sites are needed. It would also be interesting to have different levels of nitrogen input at these sites to assess more clearly the nitrogen limitation on the CO₂ fertilisation effect. We are also limited by the data the sites can provide. Due to the CO₂ fumigation over the FACE sites, daily NEE cannot be measured at the FACE sites, and therefore GPP and TER estimates cannot be derived. These variables are more directly linked to productivity than leaf area index, the variable we use in our optimisations. In future optimisations, we might also want to include more nitrogen-type variables in the cost function, such as leaf nitrogen content. There are other processes at play that also need to be assessed. For example, the effect of soil moisture and the stress response to water availability will also impact the mineralisation of organic matter and, thus, nutrient availability. Finally, structural changes do need to be made to the model to better capture the inter-annual variability of simulated NPP and LAI. This highlights how we can use FACE data to identify structural issues in models, providing an important tool for model development.

Identifying structural deficiencies is the main strength of parameter estimation – we need to be sure that we are not simulating the right model output for the wrong reasons. Indeed, if we had not been able to find a set of parameters that gave a satisfactory fit to both atmospheric regimes, this would have highlighted a missing process in the model. Ideally, we would want to calibrate under ambient conditions and test the robustness of the theory under elevated CO₂; however, given all the missing (e.g. adaptation) or highly simplified processes (nutrient limitation), using both conditions is one approach to improve the overall model behaviour while highlighting these deficiencies. Although not shown, our framework also allows us to compute the posterior parameter uncertainty, which again can be very informative for model development. We do not discuss them in this paper since our imperfect set-up (i.e. the diagonal R matrix) means the information content of the observations is overestimated in the optimisation, but we do find that the uncertainty parameters are strongly reduced in all the cases.

Although we performed a sensitivity analysis to select sensitivity model parameters, a large number of parameters were kept, some of which are indirectly impacted by the processes optimised. This can pose a risk, since changing such parameters may have an important impact elsewhere in the model and, therefore, may result in a degradation in the model–data fit against processes not considered in the optimisation. Therefore, this study is only a first step toward a more comprehensive approach with more data streams.

Nevertheless, although it is by no means exhaustive, this proof-of-concept experiment highlights the importance of manipulation experiments and the additional information they can provide for model improvement. This is the first study of a global process-based model using data in this way. With more FACE sites, these types of data could be used more consistently as part of the model optimisation procedures. Other data streams, such as the normalised difference vegetation index, solar-induced fluorescence satellite data and tree rings, could also be used to complement such optimisations, giving the best constraints on the model parameters and hence on future climate predictions.

Appendix A: Optimised parameters

All processes and equations of ORCHIDEE can be found in the different documenting publications (e.g. Krinner et al., 2005) as well as on its website (https://forge.ipsl.jussieu.fr/orchidee/wiki/Documentation/UserGuide, last access: 23 August 2023). Here, we only highlight the impacted modules, summarise the equations in which the optimised parameters feature and cite the relevant publications. Parameters optimised in the study, listed in Table 1, are coloured red in the following text.

A1 Nitrogen-related processes

The nitrogen-related parameters and their equations are thoroughly described in Vuichard et al. (2019). In this version of ORCHIDEE-CN, we prescribe leaf nitrogen concentrations. This means that the leaf $\mathrm{C}/\mathrm{N}$ ratio is fixed within a prescribed range given by two of our parameters ([CN_{leaf, min}; CN_{leaf, max}]; g_C[g_N]⁻¹). To account for the nitrogen limitation on photosynthetic activity, VC_max (the maximum rate of RuBisCO-activity-limited carboxylation) becomes a function of the leaf nitrogen content (N_l) as proposed by Kattge et al. (2009):

$$\begin{array}{}\text{(A1)}& {\text{VC}}_{\text{max}}={\text{NUE}}_{\mathrm{opt}}\cdot {\mathrm{N}}_{\mathrm{l}},\end{array}$$

where NUE_opt (µmol CO₂ s⁻¹[g${}_{{\mathrm{N}}_{\mathrm{leaf}}}$]⁻¹) is the nitrogen-use efficiency.

N_l decreases exponentially from the top to the bottom of the canopy with decreasing light levels or increasing canopy depth. The value of N_l at the top of the canopy, N_l(0), is expressed as a function of the total canopy N content, N_tot, and the LAI of the total canopy, L_tot:

$$\begin{array}{}\text{(A2)}& N_{\mathrm{l}}\left(\mathrm{0}\right)={\displaystyle \frac{K_{\text{N}}{\text{N}}_{\text{tot}}}{\mathrm{1}-\exp (-K_{\text{N}}{L}_{\text{tot}})}}.\end{array}$$

K_N is a specific extinction coefficient. Note that this is different to the extinction coefficient k used to calculate the light profile within the canopy, although both are optimised. As we decrease through the canopy, the value of N_l at a cumulative LAI (L) is defined following Dewar et al. (2012):

$$\begin{array}{}\text{(A3)}& {\mathrm{N}}_{\mathrm{l}}\left(L\right)={\mathrm{N}}_{\mathrm{l}}\left(\mathrm{0}\right)\exp (-K_{\text{N}}L).\end{array}$$

It is assumed that N_l varies through the canopy due to variations in the specific leaf area (SLA; i.e. leaf area divided by leaf mass) instead of variations in the leaf nitrogen concentration, which are kept constant. The SLA at the bottom of the canopy (SLA_init) is fixed and is also optimised.

The model calculates the nitrogen required (GN_init, g m⁻² d⁻¹) to satisfy the new carbon G_C (g m⁻² d⁻¹) to the different reservoirs under the assumption that CN_leaf does not vary (Zaehle and Friend, 2010).

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\mathrm{GN}}_{\mathrm{init}}=& ({\mathrm{FCN}}_{\mathrm{l}}/{\mathrm{CN}}_{\mathrm{leaf}}+{\mathrm{FCN}}_{\mathrm{root}}/{\mathrm{CN}}_{\mathrm{root}}\\ +& {\mathrm{FCN}}_{\text{f}}/{\mathrm{CN}}_{\mathrm{fruit}}+{\mathrm{FCN}}_{\mathrm{wood}}/{\mathrm{CN}}_{\mathrm{sap}}){\mathrm{G}}_{\mathrm{C}}\end{array}\end{array}$$

FCN_i represents the fractions (unitless) of carbon allocated to leaf (l), roots (root), fruit (f) and sapwood or stalks (wood), and CN_i represents the $\mathrm{C}/\mathrm{N}$ ratios for the different biomass pools at the previous time step. FCN_root and FCN_wood are optimised in this study. R_leaf and R_root are the fractions of N re-translocated when shedding leaves and roots (f_trans parameter in Zaehle and Friend, 2010). CTE_bact is a parameter relating denitrifier bacteria activity to soil organic matter.

Root density follows an exponential profile, with more roots located in the top soil layers. The root density profile parameter z defines the depth above which ∼65 % of roots are stored and is used to calculate plant moisture availability (Krinner et al., 2005, Eq. A18). Finally, V_maxUPTAKE is used to calculate plant N uptake (Zaehle and Friend, 2010, Supplement Eq. 8).

A2 Allocation

Allocation in ORCHIDEE-CN follows the formalisms of the OCN model (Zaehle and Friend, 2010) further described in Naudts et al. (2015) and respects the pipe model theory (Shinozaki et al., 1964). K_LAtoSA (whose range K_LAtoSA,min, K_LAtoSA,max is calibrated) is used to derive a scaling factor between leaf and sapwood mass:

$$\begin{array}{}\text{(A5)}& d_{\mathrm{l}}=K_{\text{LAtoSA}}\times m_{\text{w}}\times d_{\mathrm{s}},\end{array}$$

where d_l is the one-sided leaf area of an individual plant, d_s is the sapwood cross-sectional area of an individual plant and m_w is the water stress. Sapwood mass (M_s) and root mass (M_r) are related as follows (following Magnani et al., 2000):

$$\begin{array}{}\text{(A6)}& M_{\mathrm{s}}=K_{\text{sar}}\times d_{\mathrm{h}}\times M_{\mathrm{r}},\end{array}$$

where the parameter K_sar is calculated:

$$\begin{array}{}\text{(A7)}& K_{\text{sar}}=\sqrt{(K_{\text{root}}/K_{\text{sap}})\times (k_{\mathit{\tau }s}/k_{\mathit{\tau }r})\times k_{\mathit{\rho }s}},\end{array}$$

where K_root is the hydraulic conductivity of roots; K_sap is the hydraulic conductivity of sapwood; K_τs and K_τr are the longevity of sapwood and root, respectively; and k_ρs is the sapwood density.

A3 Phenology

For the phenology parameters, we mostly refer to MacBean et al. (2015). The photosynthetic efficiency of leaves depends on their age L_age. Using four separate age classes, biomass newly allocated to leaves goes into the first age class, and leaf biomass is then transferred from one class to the next based on a PFT-specific critical leaf age value, L_agecrit. In temperate deciduous broadleaf forests, leaf senescence is triggered when the monthly air surface temperature goes below a threshold temperature:

$$\begin{array}{}\text{(A8)}& T_{\text{threshold}}=T_{\text{senes}}+C_{\mathrm{1}}{T}_{\mathrm{l}}+C_{\mathrm{2}}{T}_{\mathrm{l}}^{\mathrm{2}},\end{array}$$

where T_l is the long-term mean annual air surface temperature and T_senes, C₁ and C₂ are PFT-dependent parameters. Once senescence has begun, a fixed turnover rate is applied, with trees losing their fine roots at the same rate as their leaves:

$$\begin{array}{}\text{(A9)}& \mathrm{\Delta }B=B\cdot \mathrm{\Delta }t/L_{\text{fall}},\end{array}$$

where Δt is the daily time step, B is the total biomass and L_fall is optimised.

A4 Photosynthesis

Stomatal conductance (g_s) is coupled to leaf photosynthesis by the following equation:

$$\begin{array}{}\text{(A10)}& g_{\mathrm{s}}=g_{\mathrm{0}}+{\displaystyle \frac{A+R_{\mathrm{d}}}{C_i-C_{i*}}}{f}_{\text{VPD}},\end{array}$$

where g₀ is the residual stomatal conductance when irradiance approaches zero, A (µmol m CO₂ m⁻² s⁻¹) is the net assimilation rate, C_i (mol CO₂ m⁻²) is the intercellular CO₂ partial pressure, C_i* is the C_i-based CO₂ compensation point in the absence of respiration (R_d) and f_VPD is the function for the approximal effect of the leaf-to-air vapour pressure difference (VPD, kPa):

$$\begin{array}{}\text{(A11)}& f_{\text{VPD}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}/(A_{\mathrm{1}}-B_{\mathrm{1}}\text{VPD})-\mathrm{1}}}.\end{array}$$

The empirical factors A₁ (unitless) and B₁ (k Pa⁻¹) are optimised in this work.

A5 Respiration

Q₁₀ (unitless) is used to calculate the temperature control of heterotrophic respiration:

$$\begin{array}{}\text{(A12)}& c_T=min(\mathrm{1},Q_{\mathrm{10}}^{(T-\mathrm{30})/\mathrm{10}}),\end{array}$$

where T is the surface or soil temperature for the above- or below-ground pools.

The growth respiration is calculated as a fraction of the remaining total biomass:

$$\begin{array}{}\text{(A13)}& R_{\mathrm{g}}={\text{FRAC}}_{\text{growthresp}}\times max(B\mathrm{-}\mathrm{\Delta }t\times \sum R_{\mathrm{m},i},\mathrm{0.2}\times B).\end{array}$$

B is the total biomass, Δt is the time step (1 d) and FRAC_growthresp is a fraction to be optimised.