Using terrestrial laser scanning to constrain forest ecosystem structure and functions in the Ecosystem Demography model (ED2.2)

Terrestrial biosphere models (TBMs) are invaluable tools for studying plant–atmosphere interactions at multiple spatial and temporal scales, as well as how global change impacts ecosystems. Yet, TBM projections suffer from large uncertainties that limit their usefulness. Forest structure drives a significant part of TBM uncertainty as it regulates key processes such as the transfer of carbon, energy, and water between the land and the atmosphere, but it remains challenging to observe and reliably represent. The poor representation of forest structure in TBMs might actually result in simulations that reproduce observed land fluxes but fail to capture carbon pools, forest composition, and demography. Recent advances in terrestrial laser scanning (TLS) offer new opportunities to capture the three-dimensional structure of the ecosystem and to transfer this information to TBMs in order to increase their accuracy. In this study, we quantified the impacts of prescribing initial conditions (tree size distribution), constraining key model parameters with observations, as well as imposing structural observations of individual trees (namely tree height, leaf area, woody biomass, and crown area) derived from TLS on the state-of-the-art Ecosystem Demography model (ED2.2) of a temperate forest site (Wytham Woods, UK). We assessed the relative contributions of initial conditions, model structure, and parameters to the overall output uncertainty by running ensemble simulations with multiple model configurations. We show that forest demography and ecosystem functions as modelled by ED2.2 are sensitive to the imposed initial state, the model parameters, and the choice of key model processes. In particular, we show that:

• Parameter uncertainty drove the overall model uncertainty, with a mean contribution of 63 % to the overall variance of simulated gross primary production.

• Model uncertainty in the gross primary production was reduced fourfold when both TLS and trait data were integrated into the model configuration.

• Land fluxes and ecosystem composition could be simultaneously and accurately simulated with physically realistic parameters when appropriate constraints were applied to critical parameters and processes.

We conclude that integrating TLS data can inform TBMs of the most adequate model structure, constrain critical parameters, and prescribe representative initial conditions. Our study also confirms the need for simultaneous observations of plant traits, structure, and state variables if we seek to improve the robustness of TBMs and reduce their overall uncertainties.

1 Introduction

Terrestrial biosphere models (TBMs) are key tools to understand the ecosystem response to anthropogenic disturbances and climate change (Medvigy and Moorcroft, 2012; McGuire et al., 2001). They are intensively used, as is or embedded in earth system models, to study plant–atmosphere interactions and predict the future of ecosystems facing global change (Poulter et al., 2010). Yet, the usefulness of TBMs is currently limited by the large uncertainties in their projections, which originate from different sources (Lin et al., 2011).

Forest structure has long been recognised as a critical component to understand forest dynamics (Hurtt et al., 2010). It influences the climatically important fluxes of carbon, energy, and water (Bonan, 2008). Yet, its realistic representation is challenging and an urgent priority in the development of next-generation TBMs (Fisher et al., 2018). The representation of the forest structure within TBMs is associated with three sources of uncertainty: model structure, model initialisation, and model parameter uncertainty.

The model structure entails, by definition, all the processes included in a model, how they are implemented, and all the underlying assumptions (Bonan, 2019). Model structure complexity varies among TBMs and also depends on the user configuration choices: different formulations of the same process can co-exist within a TBM. This complexity results from the necessary compromise between an accurate representation of reality on the one hand and the computational demand and observational requirements on the other (Shiklomanov et al., 2020). Model intercomparison studies have demonstrated that discrepancies in the representation of key processes such as forest structure (Fisher et al., 2018) or photosynthesis (Rogers et al., 2017) lead to significant uncertainties in the projections of critical variables such as the overall land carbon sequestration capacity (Friedlingstein et al., 2019, 2014, 2006; Lovenduski and Bonan, 2017).

The initialisation uncertainty reflects the error made when determining the initial conditions of the modelled ecosystem. Several approaches exist for initialising TBMs, the most common of which is probably to start runs from near-bare ground conditions, force the simulations with relevant climate forcings, and wait for the model to reach an equilibrium state: the so-called potential vegetation (Antonarakis et al., 2011). Yet, such a spin-up approach does not guarantee reliable initial demography, carbon pools, or ecosystem structure. Alternatively, forest inventories can be used to prescribe the initial composition of the ecosystem (Medvigy et al., 2009). The derivation of the initial states of critical variables, such as the aboveground biomass or the total leaf area from the plant size distribution, then relies on model default allometries, which are often derived from other, potentially non-representative, site-specific data.

Parameter uncertainty arises among other things from the necessary simplification of the natural complexity into a coherent list of model parameters, the uncertainty in the measurements used to calibrate the model, or the methods used to upscale local measurements to scales at which TBMs operate (Zaehle et al., 2005). Previous sensitivity analyses have underlined the critical importance of parameter uncertainty for the projections of ecosystem demography and productivity (Dietze et al., 2014; Massoud et al., 2019; Raczka et al., 2018; Wramneby et al., 2008). In a recent comparative study, parameter uncertainty was even shown to dominate the overall model uncertainty over process uncertainty (Shiklomanov et al., 2020). Among the model parameters, allometric coefficients scale the shape and mass of plants or their components with their size (Chave et al., 2014). Not surprisingly, multiple TBMs were shown to be sensitive to such allometric parameters (Collalti et al., 2019; Cano et al., 2020; Esprey et al., 2004). Parameter uncertainty can be reduced by constraining the range of variation of model parameters through the assimilation of different sources of observations or via model optimisation (LeBauer et al., 2013). In the past, TBMs have often been calibrated with eddy covariance data (Fer et al., 2018; Rezende et al., 2016; Collalti et al., 2016). While this approach ensures that the model correctly reproduces the short-timescale (diurnal/seasonal) dynamics of land fluxes, it does not ensure an accurate representation of forest structure and carbon pools. This is especially true because forest-structure-related parameters can present a low sensitivity to those observations (LeBauer et al., 2013; Richardson et al., 2010), and the equifinality in TBMs (Luo et al., 2009) can lead to acceptable land fluxes with a poor representation of ecosystem structure (i.e. fluxes can be reproduced from an almost infinite range of structural possibilities, some of which will be much more likely than others).

Among the different sources of observations used to reduce model uncertainties, remote sensing from various platforms (terrestrial, air- and spaceborne) has increasingly been used to monitor and understand terrestrial ecosystems (Jones and Vaughan, 2010). Lidar (light detection and ranging) data in particular have been used in the past to initialise forest biomass and constrain predictions of TBMS (Thomas et al., 2008; Hurtt et al., 2019). The recent revolution in terrestrial laser scanning (TLS, also called terrestrial lidar) provides new opportunities for constraining TBMs and reducing the uncertainties related to the vegetation structure representation (Fischer et al., 2019). The ability of TLS to measure the distance to reflecting surfaces was initially used in ecological studies to measure simple metrics such as diameter at breast height (DBH) and tree heights (Maas et al., 2008; Hopkinson et al., 2011). Since then, TLS methods have rapidly evolved to derive more complex metrics, such as the vertical profiles of the forest structure (Jupp et al., 2009; Calders et al., 2018a, b) and whole-tree volumetric assessments (Fan et al., 2020), leading to an accurate determination of forest structure across various forest types (Calders et al., 2015; Tanago et al., 2018; Takoudjou et al., 2018; Ehbrecht et al., 2017; Stiers et al., 2018; Saarinen et al., 2021). Today, the ability of TLS to accurately represent the 3D structure of forests via quantitative structure modelling (QSM, see Raumonen et al., 2013; Hackenberg et al., 2015) represents a unique opportunity to improve our understanding of forest ecosystems under changing climates (Calders et al., 2020). In particular, TLS snapshots of vegetation ecosystems could simultaneously provide important state variables to initialise TBMs and strong constraints on some critical allometric parameters, and they could help determine the most appropriate model structure for some key processes.

In this study, we evaluated the relative contributions of different sources of uncertainty (parameters, processes, initial conditions) to the overall uncertainty of multiple simulated outputs of a specific TBM, namely the Ecosystem Demography model version 2 (ED2.2). We also explored the benefits of constraining vegetation-structure-related parameters and processes using TLS for the model performance and output variability. To do so, we ran ED2.2 simulation ensembles for a temperate forest in the UK, considering different initial states for the modelled ecosystem and varying multiple model parameters and process settings with or without TLS constraints. In other words, we assessed (i) the relative importance of the model structure, initialisation, and parameter uncertainties in the ED2.2 model representation of a temperate forest and (ii) the potential added value of TLS data for vegetation modelling. To the best of our knowledge, this study is the first attempt to constrain a TBM using TLS.

2 Material and methods

2.1 Study site and data

2.1.1 Study site

Wytham Woods is a mixed deciduous forest, predominantly broadleaved, covering approximately 40 ha. It is located 5 km northwest of Oxford in southern England (Thomas et al., 2011). Owned by Oxford University, Wytham Woods has been part of the UK Environmental Change Network (ECN) and of the Smithsonian Global Earth Observatory (SIGEO) network since 1992 and 2008, respectively, and has hosted numerous ecological studies (Savill et al., 2010). The site is classified as an ancient semi-natural woodland (Hall et al., 2001), which means that the site has been continuously covered by trees through recorded history (since at least 1600), occasionally managed, and has experienced minimal intervention (i.e. no silvicultural management) since WWII (Fenn et al., 2015). Over the 1993–2008 time period, the site was characterised by a mean annual temperature of 10 ^∘C and a mean annual precipitation of 726 mm (Butt et al., 2009). The area we simulate in this study is a 1.4 ha forest plot nested within the 18 ha long-term monitoring site, part of the ForestGEO global network of forest inventory plots. This 140 m×100 m area has a local SW-coordinate (0, 100) and local NE-coordinate (140, 200) boundary. The local origin coordinate (0, 0) was located with a differential GPS at lat 51.7750579 and long −1.33904729.

Table 1 Mean values (± 1 standard deviation) of plant traits (specific leaf area (SLA) and maximum rate of carboxylation (V_c,max)) available in the TRY database for each of the five dominant species in Wytham Woods, and the local prevalences of those species (in terms of individual density and basal area). Missing traits were unavailable in TRY. The table also summarises the abundances of those five dominant species in the 1.4 ha plot in terms of absolute and relative density and basal area, as well as the PFT mapping when more than one PFT was simulated (N_PFT>1). The community-weighted means (CWM) and standard deviations (CWSD) were obtained using the basal areas as weights. Ap: Acer pseudoplatanus, Ca: Corylus avellana, Cm: Crataegus monogyna, Fe: Fraxinus excelsior, and Qr: Quercus robur.

^∗ MH: Mid-successional hardwood trees, LH: late-successional hardwood trees

Figure 1 Initial conditions in terms of tree size distribution and species composition (a), horizontal position, basal area (the size of the circles in panel (b) is proportional to the individual basal area), and species composition (b). The colour legend for species applies to both panels and is kept the same for Fig. 2 and Table 1. In the simulations, all trees were classified into single or multiple plant functional types according to the species PFT in Table 1.

2.1.2 Field inventory and terrestrial laser scanning data

The studied plot was inventoried during the summer of 2016. All trees were located, measured, and identified at the species level. The plot is largely dominated by sycamore (Acer pseudoplatanus, 65.3 % of the 815 inventoried trees in the 1.4 ha plot; see Table 1 and Figs. 1 and S1 in the Supplement), ash (Fraxinus excelsior, 10.3 % of the stems), and hazel (Corylus avellana, 8.2 % of the stems). Oaks (Quercus robur) represent a limited fraction of the woody stems (4.3 %) but disproportionately contribute (23.4 %) to the total basal area, as they mostly consist of large trees (Table 1 and Fig. 1). From the inventory, tree DBH is 24.4 cm on average (median DBH is 19.8 cm) and ranges from 2.9 to 141.2 cm.

Three-dimensional forest structure data were collected using a RIEGL VZ-400 terrestrial laser scanner (RIEGL Laser Measurement Systems GmbH) in leaf-on (June and July 2015) and leaf-off (December 2015 and January 2016) conditions (Calders et al., 2018a). The RIEGL instrument uses on-board waveform processing and records multiple return lidar data, which improves vertical sampling (Lovell et al., 2003; Calders et al., 2014). Individual trees were extracted using Treeseg (Burt et al., 2019), and their structure was modelled with TreeQSM (Raumonen et al., 2013) with the leaf-off TLS point cloud. Leaves were then added to the individual tree branches using both the leaf-off and leaf-on TLS datasets with the FaNNI algorithm (Åkerblom et al., 2018). In doing so, TLS allowed the retrieval of individual tree heights, aboveground woody biomass (modelled through estimates of volume combined with species-specific wood density), and leaf area. In addition, the individual tree crown area was computed from the vertical projection of the leaf-off point clouds of individual trees. For more details, a complete description of the TLS data collection and forest stand reconstruction is available in Calders et al. (2018a).

2.1.3 Flux tower data and species traits

Stand-scale carbon and water fluxes have been occasionally measured in Wytham Woods using the eddy covariance technique. We digitised the most recent (to our knowledge) data collected on CO₂ fluxes, which were reported by Thomas et al. (2011) for the period May 2007–April 2009. To do so, we digitised the weekly mean values of ecosystem gross primary productivity (GPP), ecosystem respiration (R_eco), and net ecosystem productivity (NEP) from Fig. 6 of the aforementioned reference using the Plot Digitizer software (v.2.6.8, http://plotdigitizer.sourceforge.net/, last access: 1 February 2021). For a more detailed description of the eddy covariance data (including the data frequency of the original data and the data quality filtering), we refer the reader to the original publication by Thomas et al. (2011).

In addition, we extracted all existing records of specific leaf area (SLA) and maximum rate of carboxylation (V_c,max) for the five most important species in Wytham Woods (Acer pseudoplatanus, Corylus avellana, Crataegus monogyna, Fraxinus excelsior, and Quercus robur) from the TRY database (Kattge et al., 2020); see Table 1 (the complete list of references from which the data originate is available in Sect. S1 in the Supplement). Individual traits were converted into ED2.2 units (m² kg${}_{\mathrm{C}}^{-\mathrm{1}}$ for SLA with a fixed leaf carbon content of 0.5 and $\mathrm{\mu }{\mathrm{mol}}_{\mathrm{C}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for V_c,max). V_c,max data were also rescaled to the ED2.2 reference temperature (15 ^∘C) using the model's default value for the temperature coefficient Q₁₀ of 2.4. Following Asner et al. (2017), we calculated the community-weighted mean (CWM) and community-weighted standard deviation (CWSD) for both traits based on the species composition and species-level average values and using the species basal areas as weights:

$$\begin{array}{}\text{(1)}& {\displaystyle }\text{CWM}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{w}_i{x}_i}{{\sum }_{i=\mathrm{1}}^N{w}_i}}\text{(2)}& {\displaystyle }\text{CWSD}=\sqrt{{\displaystyle \frac{\frac{{\sum }_{i=\mathrm{1}}^N{w}_i{\left(x_i-\text{CWM}\right)}^{\mathrm{2}}}{\left(N-\mathrm{1}\right){\sum }_{i=\mathrm{1}}^N{w}_i}}{N}}},\end{array}$$

where N is the total number of species for which data were available in TRY for each trait x, x_i is the mean trait value for species i, and w_i is the species weight (here the basal area of the species).

Flux tower data were used as a validation dataset, while the TRY data were used to constrain parameters of the TBM used in this study and described just below.

2.2 Model

2.2.1 The terrestrial biosphere model ED2.2

ED2.2 is a terrestrial biosphere model that can simulate the vegetation dynamics of a wide range of ecosystems from boreal to tropical forests (Longo et al., 2019a). It is a cohort-based, spatially implicit model that approximates the behaviour of an individual-based, spatially distributed vegetation model through a system of size- and age-structured partial differential equations (Moorcroft et al., 2001). ED2.2 integrates modules of plant growth, mortality, phenology, disturbance, hydrology, and soil biogeochemistry to predict e.g. the demography, the succession, and the dynamics of water and carbon within the simulated ecosystem.

In ED2.2, the inter- and intra-specific diversity is represented by a set of plant functional types (PFTs) that differ in leaf physiology, phenology, growth and allocation strategies, mortality, and sensitivity to environmental conditions (Medvigy et al., 2009). The trees inventoried in Wytham Woods were classified as either mid- or late-successional temperate deciduous trees (see below for the reasoning for the mapping). These PFTs are cold deciduous, i.e. leaf phenology is prognosed by the accumulation of growing degree days (growing season) and chilling days (senescing season) (Longo et al., 2019a). A comprehensive model description, including photosynthesis, allometries, radiative transfer, and phenology, is available in Longo et al. (2019a).

2.2.2 Model initialisation and forcings

In this study, the ED2.2 model was initialised using (i) near-bare ground (NBG) initial conditions (i.e. seedlings only), (ii) the field inventory, or (iii) the TLS-reconstructed size distribution. In the latter two configurations, the 1.4 ha site was initially divided into 35 square patches of 20×20 m. These three types of initial conditions are referred to below as NBG, Census, and TLS, respectively. Simulations were run for multiple years using the local forcing data for the corresponding years of the CRU-NCEP reanalysis dataset (Viovy, 2018). Simulations were run for either 5 years (Census and TLS configurations) or the approximate age since the last large-scale disturbance (100 years, NBG configuration); see Table 5. Soil texture was set according to the dominant soil type (clay) based on site-level observation (Butt et al., 2009).

2.2.3 Allometries and model parameters

In ED2.2, the carbon made available from net assimilation is partitioned at the cohort level into the different plant pools according to DBH-dependent allometries (Longo et al., 2019a). In other words, plant cohorts allocate the carbon assimilated through photosynthesis to living tissues (i.e. fine roots, sapwood, leaves, seeds), the non-structural storage pool, and the dead tissues (i.e. coarse roots, and aboveground woody biomass) depending on (i) a set of allometries and (ii) whether the plant carbon balance and environmental conditions are favourable for growth. In ED2.2, aboveground woody biomass, height, leaf biomass, and crown area are scaled through DBH-dependent allometries (Table 3). The ED2.2 default allometric models and parameters are defined according to Medvigy et al. (2009) for the leaf biomass and height, Dietze et al. (2008) for the crown area, and Albani et al. (2006) for the aboveground woody biomass.

To estimate the relative contribution of the parameter uncertainty to the variability of the model outputs, we used parameter distributions from previous ED2.2 parameter uncertainty studies (Dietze et al., 2014; Shiklomanov et al., 2020; Raczka et al., 2018; Viskari et al., 2019). We only targeted those parameters that were shown to significantly contribute to the overall parameter uncertainties in the aforementioned studies (Table 4) and set the rest to their ED2.2 default values for all simulations. For SLA and V_c,max in particular, we defined two types of parameter distributions: either relatively wide priors, as in the previous sensitivity analyses listed above (Table 4), or constrained posteriors generated by the trait meta-analysis of the Predictive Ecosystem Analyzer (PEcAn) run with the existing data in TRY and without random effects (see Raczka et al., 2018; Meunier et al., 2020; LeBauer et al., 2013). The meta-analysis was informed by TRY data only. Those distributions are referred to below as without or with TRY constraints, respectively. The uncertainty of the allometric coefficients was determined either by the range of variation of those parameters in the ED2.2 model for hardwood tree PFTs (NBG and Census configurations) or by the posterior distributions of these parameters generated when fitting the TLS data (see below).

Table 2 List of varying processes included in the model ensembles in order to evaluate the model structural uncertainty, as well as their possible configurations. Adapted from Shiklomanov et al. (2020).

Table 3 List of allometries modified in this study, ED2.2 default and TLS-derived allometric coefficients (for one or multiple simulated PFTs). The corresponding curves are plotted in Fig. 2.

^a DBH: diameter at breast height (cm); ^b MH: mid-successional hardwood trees, LH: late-successional hardwood trees.

Table 4 Description of the ED2.2 parameters varied in this study, their units, and the definitions of their priors, which were used to evaluate the model parameter uncertainty. “Source code name” is the name of the parameter as it appears in the ED2.2 source code. When trait plasticity is enabled, both SLA and V_c,max may change over time and for different cohorts of the same PFT.

^a Lnorm: log-normal distribution. ^b The values a and b define the parameters of the prior distributions (LeBauer et al., 2013). ^c NIR: near-infrared

2.2.4 Model configurations

To assess the importance of the model structure uncertainty, we targeted processes that were shown to induce significant variability in the model outputs in previous studies (Shiklomanov et al., 2020). In detail, we ran the model with multiple combinations of the following configurations: (i) closed canopies versus crowns of finite radii; (ii) two-stream versus multiple-scatter canopy radiative transfer models (RTMs); (iii) static versus plastic (varying with the available light level) SLA and V_c,max; and (iv) a single versus two plant functional types (Table 2).

By default in ED2.2, plant canopies are represented as infinitely thin flat crowns (a.k.a. complete shading or a closed canopy) that occupy virtually the entire horizontal space of the patch in which the cohort is located. In an alternative configuration, cohorts are still stacked on top of each other but have a finite radius, and hence the tallest plants only partially shade the underlying cohorts. In other words, the crown sub-model of ED2.2 determines the nature of the light competition between cohorts. Closed canopies have been shown to dramatically suppress competition from sub-dominant PFTs and typically result in unrealistically homogeneous patches (Fisher et al., 2015), while understorey cohorts receive more incoming diffuse and direct light if finite crowns are simulated.

The second sub-model we investigated was the choice of RTM. In both options (two-stream and multi-scatter), the full vertical radiation profile within each patch is resolved as a function of the canopy structure (e.g. leaf and wood area, clumping) and the environmental conditions (e.g. incident solar radiation, solar angle) following the approach of CLM 4.5 (Oleson et al., 2013). Both RTMs differ in the numerical resolution of the radiative transfers. By default (two-stream), the special multi-canopy solution of the two-stream approximation for vegetation canopies (Sellers, 1985) is used as described in Longo et al. (2019a), while the multiple-scatter is derived from first principles (Zhao and Qualls, 2005) to address the long-known issues and biases of the two-stream model (Wang, 2003). The multiple-scatter configuration increases diffuse light levels in the understorey as compared to the default two-stream approach (Shiklomanov et al., 2020).

The third sub-model that we evaluated is related to trait plasticity. By default (static), all cohorts of a given PFT share the same set of parameters, which do not evolve over time, in contradiction with the well-documented intra-specific variability of plant traits with environmental conditions (Keenan and Niinemets, 2016). In the alternative configuration (plastic), cohort SLA and V_c,max respectively decrease and increase with light availability, following empirical relationships from the tropics (Lloyd et al., 2010).

Finally, we also evaluated the impact of simulating one or multiple PFTs by either classifying all trees in the Wytham Woods inventory as belonging to the mid-successional hardwood tree PFT of ED2.2 (N_PFT=1) or according to a classification similar to the one of Dietze and Moorcroft (2011), N_PFT=2, supplemented by a clustering analysis of the allometric relationships derived from the TLS data (see below).

2.3 Analyses

2.3.1 Impact of TLS data on model allometries and initial conditions

We first compared the model default allometries with site-specific ones constrained from the TLS data. To do so, we fitted the individual plant metrics (height, crown area, aboveground woody biomass, and leaf area) versus DBH relationships derived from TLS with the set of equations used in ED2.2 (Table 2). More specifically, we fitted the parameters of the four allometries of ED2.2 using a Bayesian approach and the brms package of R (Bürkner, 2017). To account for the uncertainty of the data, we repeated the same analysis multiple times (N=100) by randomly sampling data with replacement and aggregating the resulting allometric parameter posterior distributions. To convert the leaf area obtained from TLS into leaf biomass, we used the CWM of SLA. We evaluated the quality of fit of the allometric models by computing the root-mean-square deviations or RMSD (van Breugel et al., 2011) normalised by the observed mean and the Watanabe information criterion (WAIC) for all four allometric models (height, crown area, aboveground woody biomass, leaf biomass). We fitted all allometric models using multiple possible species-to-PFT classifications and only retained the classifications that minimised the WAIC for the configurations N_PFT=1 and N_PFT>1.

To assess the relative importance of TLS for the model initialisation, we compared the tree size distributions obtained from the field inventory and the TLS data and computed the absolute and relative differences between both DBH distributions (ground truthing of TLS).

Table 5 Summary of the model configurations used in this study and the underlying model settings.

2.3.2 Ensemble runs

For each type of initial conditions (NBG, Census, or TLS), we ran ensembles of 500 simulations with parameters randomly sampled from the parameter distributions (Table 4) and with the process configuration randomly selected from the different options (Table 5). Each ensemble was equally split between runs with (250) and without (250) TRY constraints on SLA and V_c,max. The same parameter samples and process configurations were used for all three types of initial conditions and with and without TRY restrictions on SLA and V_c,max to allow independent evaluation of the impacts of the initial conditions and the TRY and TLS constraints at specific parameter values.

2.3.3 Sensitivity analyses and variance decomposition

Finally, we assessed which processes and parameters contributed the most to the overall model variance by performing a sensitivity and variance decomposition analysis following Dietze et al. (2014) and LeBauer et al. (2013). This analysis allows the fraction of the variance in target output variables attributable to individual parameters and processes (or “partial variance”) to be predicted. We chose as target output variables the ecosystem GPP during the most productive month (June) or over the leaf-on season (May–October), the total leaf area index (LAI), and the understorey photosynthetically active radiation (PAR) in leaf-on conditions, as well as the aboveground woody biomass at the end of the simulation. For the NBG configuration, we also decomposed the variance of the total stem density (which is prescribed in the other two configurations). Parameters included in the variance decomposition analyses were re-classified as belonging to one of the following three categories: allometric parameters, TRY-constrainable parameters (SLA and V_c,max), and others. All 5 years of the Census and TLS configurations were kept for analysis, while only the last 5 years of the NBG runs were considered. Note that the variance partitioning algorithm that we used only attributes to the parameters and processes their direct effects: interactions are not accounted for in the variance decomposition.

All analyses presented in this study were performed using R 3.6 (R Core Team, 2019).

3 Results

3.1 Impact of TLS data on model allometries and initial conditions

TLS-extracted and field inventory DBHs were very well correlated (R²=0.98, slope of the inventory vs. TLS linear model = 0.998; see Fig. S4). The mean (median) relative difference between the TLS and field inventory DBHs was −0.2 % (−1.7 %); see Fig. S5. The minimum and maximum absolute differences in DBH were −13.8 and 32.9 cm, respectively; the minimum and maximum relative differences were −42 % and 101 %, respectively (Fig. S5). The total tree basal area from the inventory was 36.8 cm² m⁻², while the total tree basal area obtained from TLS tree reconstruction was 36.2 cm² m⁻².

Figure 2 TLS-derived (grey, considering all tree species belonging to a single PFT) and model default (black, mid-successional hardwood trees in ED2) allometries for the aboveground woody biomass (a), tree height (b), crown area (c), and leaf biomass (d). The data to which the TLS allometries were fitted (coloured points corresponding to the tree species detailed in Fig. 1) were obtained using TLS. The coefficients used to plot the best fit and default allometries can be found in Table 3.

Individual tree measurements from QSMs applied to the TLS point cloud could all be satisfactorily represented by the ED2.2 allometric equations and a single PFT (Fig. 2). R² values of the allometric models for the individual aboveground woody biomass, height, crown area, and leaf biomass respectively reached 0.95, 0.83, 0.67, and 0.77. The normalised RMSDs changed from 18.3 % to 16.9 % (height), from 85.1 % to 75.7 % (crown area), from 146.1 % to 95.0 % (woody biomass), and from 151 % to 83.5 % (leaf biomass) when the ED2.2 default allometries for the mid-successional hardwood tree PFT were switched to TLS-derived, site-specific ones (Table 3).

Over the DBH range in Wytham Woods, TLS-derived allometries led to systematically larger allocations to aboveground woody biomass (+73 % on average; up to +177 % for the smallest tree) and leaf biomass (+75 % on average) and smaller tree heights (−1.9 m on average) as compared to the ED2.2 defaults (Fig. 2). Individual crown areas derived from TLS measurements varied between 0.2 and 465.4 m², with a mean of 26 m². As compared to the TLS-calibrated allometries, default model coefficients predicted larger crown areas for trees with DBH <64 cm (−22 % on average) and smaller crown areas for trees with DBH ≥64 cm (+17 % on average); see Fig. 2. The latter category (DBH ≥64 cm) comprised 30 trees (3.7 % of the total) and contributed 30.7 % of the total basal area and 24.9 % of the total leaf area.

Increasing the number of PFTs only slightly improved the goodness of fit of all four allometric models. The best species-to-PFT mapping according to the literature-informed minimisation of the Watanabe information criterion was to classify Acer pseudoplatanus as belonging to the late-successional hardwood PFT and the rest of the tree species as belonging to the mid-successional hardwood PFT (Table 1, Figs. S2 and S3). Using this classification, the normalised RMSDs of the allometric models decreased from 16.9 % to 16.8 % (height), from 75.7 % to 71.1 % (crown area), from 95.0 % to 77.9 % (aboveground woody biomass), and from 83.5 % to 73.9 % (leaf biomass). This mapping resulted in larger crown areas and larger carbon allocation to woody and leaf tissues for small (DBH <50 cm) trees of the mid-successional tree PFT and taller late-successional trees across all DBHs (+1.16 m on average).

3.2 Ensemble runs

Regardless of the TRY constraints and the initial conditions, the model ensembles could on average reproduce both the amplitude and the seasonality of the gross ecosystem productivity, as observed by the eddy covariance flux tower, with a maximum GPP in June and a leaf-off season with close-to-zero GPP in December–February (Fig. 3). The R² of observed vs. simulated monthly means of GPP was larger than 0.93 for all configurations (NBG, Census, TLS), while the RMSE varied between 1.2 (NBG), 1.3 (TLS), and 1.9 (Census) $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ – much lower than the mean and standard deviation of the 2 years of observational data of GPP (5.5 and 4.7 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively). Because we only simulated fully deciduous tree PFTs, model ensembles underestimated GPP during winter: simulated ecosystem LAI and hence ecosystem gross productivity dropped to almost zero in December–February (Fig. S6), while measured ecosystem productivity was non-null during the same period (Fig. 3); this was driven by evergreen understory plants such as shrubs that were not included in our simulations.

Figure 3 Seasonal cycle of the ecosystem GPP, as observed by eddy-covariance data (black dots) or as simulated by ED2.2 for multiple model configurations (columns) and with or without TRY constraints on SLA and V_c,max (rows). The thick green lines are the ensemble means, while the shaded envelopes encompass 95 % of the ensemble members. The individual ensemble members are also plotted as thin grey lines. The vertical error bars for the flux tower data represent the 95 % confidence interval of the monthly GPP. The settings of the model configurations are detailed in Table 5.

The variability of the simulated GPP was critically influenced by the model configuration and the application of constraints on SLA and V_c,max (Fig. 3). The standard deviation of the ensemble runs for the simulated GPP was, not unexpectedly, the largest for the configuration with the least information on the ecosystem (the NBG configuration without TRY constraints), and reached 6.33 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for June (Fig. 3). More than 23 % of the runs in that configuration led to unvegetated conditions (LAI <0.1 m² m${}^{-\mathrm{2},}$ all year long; see Fig. S6) after 100 years of simulations, while about 5 % of the runs simulated unrealistically dense tree cover (LAI >10 m² m⁻² in summer). Combined with the uncertainty of all other parameters, including photosynthetic ones, the LAI variability explains the extreme variability of the simulated ecosystem gross productivity. The 95 % confidence interval of the simulated ecosystem GPP in June for the NBG configuration without TRY constraints (0–19.8 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) was almost twice as large as the observed GPP at that moment (13.2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$).

Figure 4 Decomposition of the simulated GPP variance into process (orange), parameter (green), and residual (mauve) uncertainties for multiple model configurations (columns) and with or without TRY constraints on SLA and V_c,max (rows). The parameter uncertainty was further decomposed into the contributions of the allometric, TRY-constrainable (SLA and V_c,max), and other parameters (shades of green). The radii of the pie charts are proportional to the total variance of the ecosystem GPP in each configuration for the month of June (maximum GPP). The settings of the model configurations are detailed in Table 5.

Prescribing initial conditions reduced the variability of the simulated outputs: the ensemble standard deviation of GPP in June for the Census configuration without TRY constraints was 4.83 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. However, for the ecosystem productivity, constraining SLA and V_c,max was even more critical: the ensemble standard deviation of GPP in June for the Census configuration with TRY constraints decreased to 1.99 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (see Fig. 3 and also Fig. 4, where the pie chart radius is set proportional to the variance of the simulated ecosystem GPP). When both parameters were constrained and realistic initial conditions were prescribed to the model (i.e. going from the NBG without TRY constraints to the Census with TRY constraints configuration), the variability of the simulated GPP experienced a 3-fold decrease. Similarly, the variabilities of LAI (Figs. S6 and S7) and AGB (Fig. S8) were drastically reduced, with a 4-fold and a 2-fold decrease, respectively.

Given the similarities of the tree size distributions derived from the inventory and TLS (see Sect. 3.1), prescribing initial conditions had a similar impact on the variability of the outputs for the TLS and the Census configurations. Combined with the constraints on allometries, it led to a reduction in the ensemble standard deviation for GPP in June to 3.78 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the TLS configuration without TRY constraints. As for the Census configuration, constraining SLA and V_c,max with TRY data had a larger impact on the model uncertainty: the ensemble standard deviation of GPP in June for the TLS configuration with TRY constraints decreased to 1.54 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. Incrementally adding the TLS-related information to the Census with TRY constraints configuration had a positive yet more limited effect on the reduction of the model variability of GPP: the ensemble standard deviation of GPP in June was reduced by 30 % when switching from the the Census to the TLS configurations with TRY constraints. Constraining allometries with TLS had a more significant impact on LAI (Figs. S6 and S7) and AGB (Fig. S8), with a 3-fold decrease of the ensemble standard deviation from the Census with TRY constraints to the TLS with TRY constraints configurations.

All in all, the predicted variability of the ecosystem LAI and GPP was the lowest for the TLS configuration with TRY constraints: 3.79±0.50 m² m⁻² for the ensemble mean (± one standard deviation) of the ecosystem LAI (Fig. S6) and 9.86±2.89 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the ensemble mean (± one standard deviation) of the ecosystem GPP (Fig. 3), both obtained during leaf-on conditions, which compared well with independent observations (Table 6). The confidence interval of the simulated ecosystem GPP in June for the TLS configuration with TRY constraints was significantly reduced (11.8–17.6 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and much closer to the confidence interval of the observations (11.5–14.6 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). In total, the variability of the simulated GPP experienced a 4-fold decrease when parameters were constrained, realistic initial conditions were prescribed, and TLS data were used to constrain the allometries (i.e. going from the NBG without TRY constraints to the TLS with TRY constraints configuration).

Table 6 Summary of most important states and fluxes in all three model configurations and how they compare with observational datasets, including flux tower data on ecosystem respiration and net ecosystem productivity. Those numbers take into account the full 5 years of simulation for the prescribed model configurations (Census and TLS), the last 5 years of simulation for near-bare-ground conditions (NBG), and the 2 years of eddy covariance observational data. For the observations of LAI in the leaf-on season, we provide a range of variation. LAI: leaf area index, AGB: aboveground biomass, GPP: gross primary production, NEP: net ecosystem productivity, PAR: photosynthetically active radiation.

^a References: Thomas et al. (2011) and Fenn et al. (2015). ^b Reference: Roberts et al. (1999).

3.3 Variance decomposition and sensitivity analyses

The variance of the ecosystem GPP was dominantly driven by the parameter uncertainty, regardless of the configuration and the application of TRY constraints (Fig. 4). Together, TRY-constrainable parameters, allometric coefficients, and the other ED2.2 parameters included in the sensitivity analysis contributed on average 63 % of the total variance of GPP in June. Constraining SLA and V_c,max with TRY datasets dramatically decreased the relative contribution of these two parameters to the overall variance: moving from uninformed priors to posteriors generated by the trait meta-analysis of PEcAn made the sum of their partial variances drop from a majority (57 % on average for all three configurations) to a small contribution (7 % on average for all three configurations), their share being mainly replaced by unconstrained parameters, which increased from 6 % to 50 % on average across all configurations (Fig. 4), especially the parameters quant. eff., clumping, and growth (Fig. 5). The variance decomposition of the simulated ecosystem LAI and aboveground biomass led to very similar results but with a larger contribution from allometric parameters: allometric parameters contributed on average 6 % and 20 % of the variances of LAI and AGB, respectively, which are larger than their contribution to the variance of GPP (3 %), illustrating the importance of TLS for constraining the ecosystem structure (Figs. 5, S7 and S8).

Figure 5 Contributions of individual or allometric parameters (Bd, Bl, CA, and height include all parameters for the respective allometries; see Table 2) to the predicted uncertainties in ED2.2 of multiple state variables (PAR: photosynthetically active radiation reaching the ground, LAI: leaf-on ecosystem leaf area index, AGB: final ecosystem aboveground biomass, GPP: leaf-on ecosystem gross primary production) for multiple model configurations (columns) with or without TRY constraints on SLA and V_c,max (rows). Only those parameters that contributed at least once to 5 % or more of the total variance were included in the panels. Parameter descriptions and distributions are given in Table 4. The settings of the model configurations are detailed in Table 5.

On average, processes only accounted for 12 % of the overall variance of GPP with the maximum (minimum) obtained for the TLS configuration with TRY constraints (the NBG without TRY constraints): 20 % (5 %). Process uncertainty was dominated by the type of crown model (5 %) and the radiative transfer model (4 %). Trait plasticity only contributed marginally to the overall variance (<1 % on average). Processes (especially the choice of the RTM) played a stronger role in the available light in the understorey (on average 40 % of the total variance), especially in runs with prescribed initial conditions (on average 56 % of the total variance; see Fig. S9). Due to compensatory effects (Fig. S2), the number of simulated PFTs had a limited impact on all of the considered model outputs: N_PFT only contributed 3 % of the variance of ecosystem GPP, 2 % of the variance of LAI and PAR, and 1 % of the variance of AGB.

Figure 6 Tree size distributions for multiple model configurations starting from near bare-ground conditions after 100 years of simulations (coloured bars), and how they compare to the field inventory (grey). The histograms and the vertical error bars represent the mean ± one standard deviation of the ensemble member runs. Only runs that generated vegetation were kept when plotting this figure.

Figure 7 Ecosystem average of the leaf area density vertical distribution for the month of June for different model configurations (coloured lines and envelopes) without (left) and with (right) TRY constraints on SLA and V_c,max. The envelopes encompass the mean ± one standard deviation of the ensemble member runs. Only runs that generated vegetation were kept when plotting the NBG envelopes. The settings of the model configurations are detailed in Table 5.

3.4 Ecosystem structure and functions

Despite similar seasonal cycles of ecosystem productivity (Fig. 3), ensemble means exhibited highly contrasting ecosystem structures (Figs. 6 and 7). None of the unprescribed simulations (NBG configuration) could capture the size distribution observed through the inventory (Fig. 6). Small-size stem (especially DBH <50 cm) densities were underestimated while large tree (DBH >100 cm) densities were overestimated in the vegetated simulations (LAI >0.1 m² m⁻²) of the NBG configuration with or without TRY constraints. Switching from closed canopy to finite crowns systematically increased the density of small (DBH <50 cm) trees by 73 % on average, just like constraining SLA and V_c,max with TRY data. While the ecosystem LAI of the NBG configuration with closed canopies compared well with independent observations from the literature (3.83±1.94 versus the range 3.6–4.1 m² m⁻² observed in Wytham Woods; Table 6), the vertical arrangement of the leaves significantly differed from what was observed by TLS and imposed in the TLS configuration (Fig. 7) as a result of the differences in tree size distribution (Fig. 6).

Despite lower total leaf areas, the infinitely wide crown configuration (closed canopies, Table 6) made the forest more opaque to the incoming solar radiation than the finite crowns. Across all configurations, the PAR available in the understory decreased by 15 % throughout the year while the ecosystem LAI decreased by 18 % when closed canopies were simulated (Table 6). For near-bare-ground configurations, the LAI of the potential vegetation simulated was 23 % lower with infinite crowns, and 16 % less PAR reached the understorey.

As the soil received more radiation when finite crowns were simulated, it was warmer and, as a result, heterotrophic (and ecosystem; see Table 6) respiration increased (+25 % on average) when switching from infinite to finite crowns. Forest carbon stocks also diverged between configurations: driven by higher allocations to leaf and aboveground woody biomass (Fig. 2), aboveground carbon storage was larger (+74 % on average) in TLS-derived runs than when default allometries were applied (Table 6). Aboveground woody biomass from configurations starting from near-bare-ground conditions was systematically underestimated compared to the TLS estimates (11.4 kg_C m⁻² on average for the NBG configuration versus 24.5 kg_C m⁻² on average for the TLS configuration). However, the larger allocation to woody biomass induced by the use of TLS-derived allometries mostly did not impact any other model outputs (Fig. 5), as that carbon pool is inert and does not influence a lot of processes downstream (e.g. more woody biomass does not translate into exacerbated light interception). Leaf biomass allometry derived from TLS both reduced the simulated LAI and ecosystem GPP to more realistic values and constrained its variability (Figs. 3 and S6 and Table 6).

None of the simulations/configurations could accurately represent all the features of Wytham Woods. The model simulations starting from near-bare-ground conditions failed to capture the vertical distribution of leaves (Fig. 6) and the tree size distribution (Fig. 7); the model simulations prescribed with the inventory overestimated the ecosystem GPP (Table 6); and the model simulations from the three configurations all overestimated the net ecosystem productivity (NEP), due to an overestimation of GPP (Census) and/or an underestimation of the ecosystem respiration (Census, NBG, and TLS); see Table 6. Model simulations underestimated R_eco on average by −17 %, leading to unrealistic NEP predictions, which illustrates the need to constrain or optimise autotrophic and heterotrophic respiration parameters along with the photosynthetic and allometric parameters to align them with observational data.

4 Discussion

4.1 The relative weights of the different sources of uncertainty

The different model configurations tested in this study led to contrasting predictions of vegetation states. Depending on the chosen model outputs, the relative weights of the sources of uncertainty considerably varied. Near-bare-ground simulations generated potential vegetations that significantly differed in their demography from observations (Fig. 4), while prescribing the initial tree size distribution was not a guarantee of accurate reproduction of observed land fluxes (Fig. 3, Table 6). The finite crown area representation also had a substantial impact on the model outputs. In particular, limiting the crown radius to finite values promoted smaller plants in the understorey (Fig. 6), increased the simulated LAI (Table 6), and profoundly modified the vertical distribution of light in the canopy (Fig. 8 and Table 6). Carbon pools also considerably diverged between model configurations, especially when TLS-derived allometries were taken into account (Table 6).

However, in general, it was the parameter uncertainty that dominated the overall model uncertainty (Figs. 3, S7 and S8), as previously also observed for ED2.2 simulations of temperate forests (Shiklomanov et al., 2020). The parameters that dominated the variance depended on the use of TRY and/or TLS constraints. When observations were available, uncertainty was transferred to other unconstrained parameters while the overall variance was reduced, as in similar previous studies (Meunier et al., 2020), which supports the need to progressively integrate observations of the most sensitive parameters until the model variance is reduced to satisfactory levels in an efficient data-model fusion loop (Dietze et al., 2014).

Although the parameter uncertainty was larger in magnitude than the process uncertainty, the crown size representation and the choice of RTMs appear to drive a significant part of the model process uncertainty and should receive more attention in future analyses, particularly because the implementation and the sensitivity of the radiative transfer processes are currently overlooked in ED2.2, like other vegetation models (Fisher et al., 2018; Viskari et al., 2019).

4.2 The added value of TLS for vegetation modelling

The quantitative information that remote sensing generates at unprecedented spatial and temporal scales can serve the purpose of reducing uncertainties in TBM projections. It has already been shown that airborne laser scanning (ALS) combined with an individual-based forest model could offer new insights into the contribution of plant size to ecosystem functioning (Fischer et al., 2019). Similarly, ALS and synthetic-aperture radar have successfully been applied to prescribe the initial structure and composition of tropical forests (Antonarakis et al., 2011, 2014; Longo et al., 2020), and lidar data have been coupled to allometric models to estimate carbon stocks and fluxes at a large scale (Hurtt et al., 2019; Thomas et al., 2008). Yet, our study is the first attempt to inform a TBM with TLS data. As compared to ALS, TLS offers a few significant advantages, as well as some drawbacks that are important to remember. Airborne techniques allow for wall-to-wall coverage characterising the 3D forest structure at the regional scale, whereas TLS offers far more detailed information but only at the local (up to a few ha) scale. However, TLS is capable of estimating the volume of individual trees directly, instead of relying on allometries that require calibration and thus field measurements. In addition, it can accurately capture the entire size distribution (DBH and height) of the sample plot, while smaller trees can easily be missed with airborne surveys (Wang et al., 2016), leading to incorrect demography, especially in dense forests.

Because TLS data are complementary to the datasets that are frequently used for model calibration (e.g. eddy covariance data), they can contribute in a collective effort towards realistic representations of ecosystems in TBMs. TLS has the potential to fill important parameter and process gaps and, in doing so, to help reduce the uncertainties in vegetation model simulations. The steep increase in the amount of available forest TLS data over the past decade (Calders et al., 2020) makes its coupling with TBMs even more timely. As demonstrated in this study, TLS observation can ensure a more adequate model structure, constrain model allometric parameters, and prescribe representative initial conditions. Yet, only a combination of constraints on both allometries (using TLS data) and photosynthetic parameters (thanks to TRY data) could satisfactorily reduce the model uncertainties to their lowest levels, which supports the integration of multiple data sources into TBMs for more realistic simulations (Peylin et al., 2016). Such a combination of a TBM and multiple data streams allowed us to accurately simulate both ecosystem productivity and ecosystem community composition with physically realistic parameters, which was previously highlighted as a challenge for dynamic vegetation models (Shiklomanov et al., 2020; Fisher et al., 2010).

In the future, TLS could inform vegetation models even more. The TLS community is indeed actively working on the derivation of additional tree- or stand-scale parameters from lidar raw data and 3D point clouds. Those parameters include leaf angle distributions (Vicari et al., 2019), clumping (Zhao et al., 2012), and reflectance (Calders et al., 2017), which have been shown to significantly contribute to the overall model uncertainty (Meunier et al., 2022; Shiklomanov et al., 2020; Viskari et al., 2019). Yet, theoretical, technological, and technical challenges specific to each parameter still need to be overcome before these sensitive traits can be constrained with TLS in a study similar to this one.

4.3 Model equifinality

Some runs from all three configurations (prescribed, or not, with the initial size distributions) could reproduce the seasonal cycle of GPP observed by the flux tower (Fig. 3). However, those “optimal” simulations were very different from a forest structure point of view (Table 6, Figs. 6 and 7). This situation illustrates the low identifiability of numerous TBM parameters and the need for multiple simultaneous constraints and observations. While knowledge of the aboveground carbon storage is critical to estimate forest sink strength and the overall carbon storage capacity of the ecosystem (Keeling and Phillips, 2007), it has a limited impact on simulated land fluxes (GPP in particular; see Fig. 5) that are often used to calibrate TBMs. The parameters controlling land fluxes, namely those controlling the ecosystem LAI (Williams and Torn, 2015; Wei et al., 2013) and those related to photosynthesis (Fig. 5), are also confounded, echoing observed trade-offs of the leaf economic spectrum (Wright et al., 2004; Peaucelle et al., 2019). TLS has the potential to discriminate equifinal model simulations with similar land fluxes but contrasting structures. On-site trait measurements (Fig. 3) could further help avoid those risks of equifinality (Babst et al., 2020; Peaucelle et al., 2019).

4.4 Study limitations

Our findings come with several important limitations. First, the eddy covariance flux data (2007–2009) preceded the observation of the forest structure (the TLS and field inventory occurred over the 2015–2016 period) by almost a decade. The forest composition and demography might have changed in the meanwhile, which reduces the confidence in the validation with eddy covariance data (Fig. 3). This is even more true as one realises that the validation dataset is rather limited in size and information content (very low year-to-year variability in observed fluxes). Yet, in this study we were more interested in the variance decomposition for different model configurations (Figs. 3 and 4) than the actual goodness of fit of every single configuration. In addition, in the absence of locally observed meteorological drivers, we had to force the model simulations with regional datasets that cannot serve the purpose of capturing the day-to-day variability or the diel cycle, which forced us to only compare the modelled and observed seasonal GPP cycle. Furthermore, GPP is not a directly observed but rather a derived (modelled) quantity, as opposed to the net ecosystem exchange of carbon and the latent heat flux of water, which are directly measured. We could not access water flux raw data; nor were they reported in publications that we knew of. GPP uncertainties were also not quantified in the original publication of Thomas et al. (2011). While NEP values were reported, validating the model simulations with those values would have biassed our analyses, as we could not constrain respiration parameters with data. Mismatches between different data sources and/or the low availability of good-quality data are recurrent issues in vegetation modelling exercises. Despite multiple initiatives to standardise high-quality data, such as Fluxnet (Baldocchi et al., 2001), we emphasise here the need for concomitant observations in experimental and observational plots.

Second, the comparison between the potential vegetations as simulated by ED2.2 and the field inventory data is also imperfect, as Wytham Woods is a managed forest that has been frequently coppiced and pollarded. The disturbance history experienced by the ecosystem is mostly unknown, preventing us from reproducing the current forest demography by the model.

Third, the trait meta-analysis was run with random effects turned off, which can generate parameter posterior distributions that are too narrow (Raczka et al., 2018) and hence underestimate the contribution of the TRY-constrained parameters (see e.g. Fig. 4). A similar analysis that includes random effects should be performed to evaluate such an underestimation.

Finally, the ecosystem growth form complexity was neglected in this study. We only simulated tree PFTs, while shrubs and grass species also coexist in Wytham Woods. Integrating this ecological complexity would not have brought additional information or robustness regarding the objectives of our study on the variance decomposition, while it would increase the dimensionality and complexity of the problem. Future research should investigate whether the main findings highlighted in this study hold with other PFTs across other sites and biomes, or even in other vegetation models (Dokoohaki et al., 2021).

5 Conclusion

Vegetation models are important tools to predict the fates of ecosystems in a changing climate, but are often used as black-box tools due to their complexity. They have been designed to realistically represent the ecosystem that they simulate but often fail to do so, primarily because of considerable parameter uncertainties as well as process and initialisation errors. Even for the state-of-the-art process-based terrestrial biosphere models, not all parameters can be constrained with data: some cannot be observed in the field or require calibration, or the appropriate observational trait data may be missing. In addition, model initialisation and the choice of model structure necessarily lead to additional uncertainties. We demonstrate in this study that TLS has the potential to provide initial condition estimates and to constrain some critical vegetation model parameters (allometries) and processes (crown representation). Combined with trait-based constraints on a few key parameters, TLS was able to define a model configuration that could reproduce both the ecosystem productivity and the plant community composition of the simulated site with physically realistic parameters, as well as to considerably reduce model uncertainties.