PDG-Arena: an ecophysiological model for characterizing tree-tree interactions in heterogeneous and mixed stands

Abstract

Key message

PDG-Arena, a new individual- and process-based ecophysiological model, was developed to simulate the functioning of mixed-species forests. It was evaluated using annual growth data on beech-fir stands of the French Prealps.

Context

In the context of the ongoing climate and biodiversity crises, mixed forest stands are increasingly considered as a sustainable management alternative to monospecific and even-aged stands.

Aims

We developed a new individual- and process-based forest growth model, PDG-Arena, to simulate mixed forest growth and functioning, and test ecophysiological interactions among trees in mixed stands.

Methods

The model builds upon the validated ecophysiological stand-scale model CASTANEA and integrates tree competition for light and water. We evaluated the performance of PDG-Arena by comparing the simulated growth with annual radial growth data from 37 common beech and silver fir monospecific and mixed plots in the French Prealps.

Results

PDG-Arena performed slightly better than CASTANEA when simulating even-age and monospecific forests (r² of 32.1 versus 29.5%). When using structure-diverse and species-diverse inventories, PDG-Arena performed better than CASTANEA in pure beech (38.3 versus 22.9%) and mixed stands (40.5 versus 36.3%), but not in pure fir stands (39.8 versus 42.0%). The new model also showed a significant positive effect of species mixing on gross primary production (+ 5.5%), canopy absorbance (+ 11.1%), and transpiration (+ 15.8%) in the tested stands.

Conclusions

Our results show that tree-level process-based models such as PDG-Arena, formally simulating interspecific interactions, can serve as a valuable tool to understand and simulate the carbon, radiative, and water dynamics of mixed stands.

1 Introduction

Understanding how forest ecosystems function is a crucial step for developing forest management strategies adapted to the challenges of climate change (Bonan, 2008; Lindner et al., 2010; Trumbore et al., 2015) and more generally global change (González de Andrés, 2019). In this context, mixed forests, in comparison with monospecific stands, have received increasing attention due to their documented ability to maintain key ecosystem services while enhancing stand resilience (van der Plas, 2016; Seynave et al., 2018; Messier et al., 2022; del Río et al., 2022).

However, the ecophysiological functioning of mixed stands is still poorly understood (Forrester, 2014; Forrester and Bauhus, 2016). In particular, even though species mixing seems on average to increase stand productivity in comparison to monospecific stands (an observation known as overyielding) (Liang et al., 2016; Zhang et al., 2012; Vilà et al., 2007; Forrester and Bauhus, 2016; Piotto, 2008), this trend depends on stand structure and species composition (Zhang et al., 2012; Ratcliffe et al., 2015), as well as abiotic conditions (Ratcliffe et al., 2016; Toïgo et al., 2015). Regarding the effect of tree species richness on the resistance of stands to drought episodes, the literature shows heterogeneous results (Grossiord, 2019). Indeed, the direction of the effect seems to depend on mixture composition—and particularly on the respective species strategies in reaction to soil water deficit (Pretzsch et al., 2013; Mas et al., 2024; Jourdan et al., 2020)—as well as on environmental conditions (Grossiord et al., 2014; Forrester et al., 2016; Pardos et al., 2021).

Stand structure, particularly tree density and size variability, can act as a confounding factor in the diversity-functioning relationship (Metz et al., 2016; Dănescu et al., 2016; Cordonnier et al., 2019; Zeller and Pretzsch, 2019). To better understand the processes underlying these relationships, it is therefore important to separate the effects of mixing related to differences in stand structure (age, size, diameter) from those related to differences in the physiology of species (crown architecture, water and nutrient use, etc.; see Forrester and Bauhus, 2016).

Furthermore, the interactions observed in a mixture may be of various kinds (Forrester et al., 2016), which could give rise to contradictory effects. For example, an increase in the amount of irradiance intercepted by tree canopies in mixtures—e.g., through crown complementarity and plasticity, see Jucker et al. (2015)—could lead to an increase in gross primary production, but also in transpiration, with a potentially negative effect on available soil water (Jucker et al., 2014). Forrester (2014) proposed a conceptual model to account for the mechanisms of interaction between diversity, functioning, and environment. In this framework, interspecific interactions resulting in reduced competition for a given type of resource generate beneficial effects for individuals when this resource becomes scarce.

Assessing and predicting the functioning of mixed stands therefore requires detailed knowledge of interspecific interactions. This knowledge must be based on interactions between individuals and on the ecophysiological processes underlying these interactions, i.e., competition for light, water, and nutrients (Pretzsch et al., 2017; Grossiord, 2019). This knowledge is all the more needed as abiotic and biotic conditions are changing due to global change (Ammer, 2019).

Although experimental and observational systems are necessary for studying the diversity-functioning relationship in forests, they are limited by sample size, number of measured variables and number of confounding factors that can be controlled (Bauhus et al., 2017). Modeling can virtually overcome these limitations, subject to the assumptions contained in the model, which depend to a large extent on our ecological knowledge as well as on the availability of climatic, pedological, silvicultural and physiological data. The modeling approach has been used to put forward hypotheses to explain overyielding in mixed forests. For example, Morin et al. (2011) showed with simulations that it could be caused by the diversity of species traits related to shade-tolerance, maximum height, and growth rate (although other explanations could not be ruled out). Simulations also make it possible to virtually assess the stability of the productivity of forest mixtures while testing numerous community compositions (Morin et al., 2014), even under unprecedented climatic conditions (Jourdan et al., 2021).

The literature (Korzukhin et al., 1996; Cuddington et al., 2013; Morin et al., 2021) depicts a spectrum ranging from empirical models, which are based on relationships calibrated from observations between final variables such as productivity and explanatory variables (e.g., rainfall, irradiance), to process-based models where final variables are computed using explicit elementary processes (e.g., photosynthesis, transpiration, phenology). For some authors (Fontes et al., 2010; Cuddington et al., 2013; Korzukhin et al., 1996), process-based models seem more relevant for simulating ecosystem functioning undergoing climate change because they can theoretically be applied to a larger range of environmental conditions than empirical ones. As a result, they now play an important role in research on the ecophysiological functioning and prediction of forest dynamics (Gonçalves et al., 2021; Barbosa et al., 2023). However, compared to empirical models, process-based models are more difficult to parameterize and rely on more assumptions about the ecological functioning of forests (e.g., the hypothesis that growth is primarily driven by photosynthetic activity, Fatichi et al., 2014). When it comes to simulating mixed stands, models that simulate elementary processes are expected to reproduce the mechanisms that lead to interspecific interactions, bringing us closer to understanding them (Forrester and Bauhus, 2016).

Among process-based models, a distinction is made between individual-based models, e.g., Jonard et al. (2020), and stand-scale models, e.g., Dufrêne et al. (2005). Several diversity-functioning studies in forests have highlighted the importance of tree-tree interactions in defining the nature of interspecific interactions at stand level (Trogisch et al., 2021; Jourdan et al., 2020; Guillemot et al., 2020; Jucker et al., 2015). Thus, the individual tree scale seems relevant for representing the key mechanisms that govern the functioning of mixed forests (Porté and Bartelink, 2002). Finally, process- and individual-based models have the ability to distinguish the effects of competition between individuals of different species (mixing effect) and the effects of competition between individuals of different sizes (structure effect). So far, few models are able to simulate mixed stands by taking advantage of both physiological mechanisms and the individual scale (Reyer, 2015; Pretzsch et al., 2015).

Here we present a new individual- and process-based forest growth model, PDG-Arena (the arena represents the stand, a place where trees compete and more generally interact). Our model was developed to observe the stand scale properties that emerge when trees of different species and size compete in a given environment. It was therefore built: (i) from elementary physiological processes using the stand-scale model CASTANEA (Dufrêne et al., 2005) and (ii) by integrating interactions among trees, notably competition for light and water.

The performance of PDG-Arena was evaluated using annual growth data from a monitoring network of monospecific and multispecific stands of common beech (Fagus sylvatica L.) and silver fir (Abies alba Mill.). Firstly, we tested whether PDG-Arena, despite increased complexity, accurately reproduces the performance of CASTANEA when both models are run under comparable conditions. Secondly, we evaluated PDG-Arena’s performance in different conditions in terms of stand structure and species diversity. Lastly, using PDG-Arena, we evaluated the effect of species mixing on carbon assimilation, water use and irradiance interception.

2 Materials and methods

2.1 Model description

2.1.1 From CASTANEA to PDG-Arena

PDG-Arena was designed as an extension of PDG (which stands for Physio-Demo-Genetics, Oddou-Muratorio and Davi, 2014), an individual-based and spatially explicit model that combines: (1) the process-based model CASTANEA to simulate tree ecophysiology, (2) demographic processes allowing to model tree survival and reproduction and (3) a quantitative genetics simulation module accounting for the heritability and intraspecific diversity of key life history traits of the CASTANEA model. While PDG is built with the idea of simulating the evolutionary dynamics of functional traits of importance for adaptive forestry in regular monospecific stands (Lefèvre et al., 2014), PDG-Arena is designed to simulate ecological interactions between trees in diverse, multispecific stands.

CASTANEA is an ecophysiological forest growth model that simulates the dynamics of homogeneous stands (Fig. 1a). Among others, it has been parameterized and validated on common beech (Fagus sylvatica L., Dufrêne et al., 2005) and silver fir (Abies alba Mill., Davi and Cailleret, 2017). CASTANEA is composed of five equal-sized leaf layers that perform photosynthesis based on stomatal conductance and on the level of radiation received by each layer, which is determined using a horizontally homogeneous multi-layer radiation model. The resulting gross primary production, minus autotrophic respiration, is then allocated into the leaf, fine root, coarse root, branch, trunk and reserves compartments (Davi et al., 2009). The amount of leaf transpiration is controlled by net radiation, stomatal conductance as well as ambient temperature and vapor pressure deficit. Stomatal conductance, limiting photosynthesis and transpiration, is controlled by soil water deficit (using the critical threshold of relative extractable water $REW_c$; Granier et al., 1999). Lastly, leaf surface growth is controlled by day length and mean temperature. The temporal scale of the processes in CASTANEA is the same as that of PDG-Arena, as shown in Table 1.

Fig. 1

Conceptual diagram of the (a) CASTANEA and (b) PDG-Arena forest growth models input and functioning. CASTANEA and PDG-Arena respectively simulate the growth of regular monospecific stands and (potentially) diverse multispecific stands. In CASTANEA, all processes, including radiation balance, carbon fluxes, transpiration of trees and soil water budget occur at the stand level, on horizontally homogeneous leaf layers. PDG-Arena takes advantage of CASTANEA carbon and transpiration processes but performs them at the tree level, while a water budget is computed at the stand level. Its radiative balance is handled by the SamsaraLight library which casts sun rays through a 3D representation of tree crowns. Processes involving competition between trees in PDG-Arena are shown in dashed boxes

Table 1 Temporal and spatial scales of physical and physiological processes in PDG-Arena

The existing model PDG considers isolated abstract trees, simulating the dynamics of each of them using stand-scale processes of CASTANEA. All quantitative physiological variables in CASTANEA and in PDG are expressed on a per area of soil basis: e.g., the gross primary production is expressed in gC m⁻². The first improvement of PDG-Arena over PDG is that the physiological processes are simulated at tree and not at stand level (Fig. 1b). To do so, physiological processes are related to the projected area of the individual crowns rather than to the stand area. This paradigm shift implied changing the definition of some variables. As depicted in Fig. 2, Leaf Area Index (LAI) is now defined for each tree as the amount of leaf surface of a tree per m² of soil below its crown. While stand LAI in CASTANEA depends on gap fraction, individual tree LAI in PDG-Arena does not: the LAI of a tree only accounts for its leaf surface and its crown projection area. The same reasoning applies to other physiological variables, such as carbon uptake, transpiration, intercepted irradiance, etc. Also, the Leaf Mass per Area (LMA), as it depends on the amount of irradiance intercepted by neighboring trees, is computed at the individual level in PDG-Arena according to the vertical profile of the leaf area of neighboring trees (see the “Computing leaf mass per area” section in Appendix 2).

Fig. 2

Representation of the Leaf Area Index (LAI) in PDG-Arena. In the case of the stand-scale model CASTANEA, LAI does consider the uncovered soil surface (a). In the individual-based model PDG-Arena, only the soil surface under the crown of the individual is considered. Values of leaf surface, soil surface, and LAI are arbitrary

The second improvement of PDG-Arena over PDG is that it integrates mechanisms of competition for light and water between neighboring trees (see Fig. 1b) by (i) making trees share the same stand soil water pool and (ii) simulating irradiance at tree level using a ray tracing model.

2.1.2 Competition for water

Competition for water is a crucial element in the dynamics of mixed stands. We modeled competition for water symmetrically between individuals, i.e., trees in the same plot all draw from the same water reservoir without spatial differentiation, either horizontal (distance between individuals) or vertical (depth).

The precipitation intensity is determined at the daily level from the climate file. A fraction of total precipitation (given by the proportion of soil that is directly under any tree crow) is distributed among trees according to their respective total leaf area. For each tree, a calculation of rainfall interception, runoff along the trunks and transmitted rainfall is performed, using the same equation than CASTANEA (Dufrêne et al., 2005). The fraction of total precipitation that was not distributed to the trees falls directly to the ground. Transpiration and crown evaporation (referred as crown evapotranspiration) are calculated individually at hourly time steps using the Penman-Monteith equation (Monteith, 1965), and taking into account the energy absorbed by individual crowns (see Sect. 2.1.3). Additionally, soil evaporation is computed hourly and homogeneously along the plot, following equations of CASTANEA (Dufrêne et al., 2005) and using the energy absorbed by the soil (see Sect. 2.1.3). Evapotranspiration from understorey vegetation is ommited.

Considering direct precipitation, runoff along trunks and transmitted rainfall on the one hand, and tree evapotranspiration, soil evaporation and drainage on the other, a water balance is computed each day at the stand level (Table 1 and Fig. 1b). Therefore, soil water status (soil moisture, litter moisture and soil potential) is the same for every tree within a plot on any given day.

2.1.3 Competition for light

Competition for light in PDG-Arena is simulated using SamsaraLight, a ray tracing library derived from Courbaud et al. (2003) and maintained on the Capsis modeling platform. The integration of SamsaraLight with the physiological model CASTANEA (which is partly inspired from the approach in the HETEROFOR model, Jonard et al., 2020) is described here. The stand soil is subdivided into square cells of 1.5 m width. Irradiance is evaluated both in the PAR (photosynthetically active radiation) and in the NIR (near infrared radiation) domains. For each domain, SamsaraLight generates a set of diffuse and direct beams, and computes their interception by either tree crowns or soil cells. The simulated energy absorbed by crowns and cells (in MJ year⁻¹) is then temporally distributed at the hourly scale. Finally, the energy absorbed by a crown is distributed among its five leaf layers, which are part of the CASTANEA model for each tree.

Definition of crowns

Each tree is represented by a crown occupying a volume in space and is defined by the following variables:

• The height of the tree h;

• Its crown base height, hcb;

• Its crown radius crownRadius;

• Its shape, which can be conical, cylindrical or ellipsoidal (a tree shape is vertically bounded by h and hcb and horizontally bounded by crownRadius);

• Its leaf area density at full vegetation, denoted LAD, in m² of leaf per m³ of crown volume;

• Its attenuation coefficient k;

• Its clumping index $\Omega$ defining the aggregation of leaves inside the crown.

Crown shapes are conical in the case of silver fir and ellipsoidal in the case of common beech. Tree h and hcb values are model inputs (see Sect. 2.2). Tree crown radii have been estimated using an allometric relationship based on species and diameter at breast height (DBH):

$$\begin{aligned} crownRadius = \beta _{crown} + \alpha _{crown} \times DBH \end{aligned}$$

(1)

$\alpha _{crown}$ and $\beta _{crown}$ are species dependent parameters estimated on site at Mont Ventoux (unpublished data from one of the authors). $\Omega$ is species dependent and was measured on Mont Ventoux sites by Davi et al. (2008). The attenuation coefficient k depends on species, radiation domain, type of radiation (direct, diffuse) and beam height angle. Its value is determined using reverse-engineering of SAIL (the radiation sub-model in CASTANEA) as described in Sect. 2.

The LAD of a tree is the ratio of its leaf area to its crown volume. The leaf area of a given tree i (denoted $LA_i$) is determined using the stand leaf area at full vegetation ($LA_{stand}$, which is a simulation input, see Sect. 2.2). For every tree, its fraction of leaf area over stand leaf area is proportional to its theoretical leaf area $LA_{th}$:

$$\begin{aligned} LA_i = LA_{stand} \times \frac{ LA_{th}(DBH_i, species_i) }{\sum _j^{n} LA_{th}(DBH_j, species_j) } \end{aligned}$$

(2)

$LA_{th}$ is given by an allometric equation based on DBH and species from Forrester et al. (2017), which covers most parameterized species of CASTANEA:

$$\begin{aligned} LA_{th}(DBH_i, species_i) = \beta _0(species_i) \times DBH^{\beta _1(species_i)} \end{aligned}$$

(3)

During the radiation balance computation, each tree LAD is at its maximum. However, a fraction of the absorbed radiations per tree is removed daily depending on their current phenological state (see the “Reduction of absorbed radiations in SamsaraLight” section in Appendix 2).

Ray casting

SamsaraLight generates two sets of beams. Firstly, diffuse rays are generated in all directions, using a 5° discretization. Secondly, direct rays are generated to follow the hourly trajectory of the sun for one virtual day per month. Each set of beams contains the energy of the entire year (in MJ) for both diffuse and direct radiations. The stand plot is subdivided into square cells of 1.5 m width. All beams are replicated for each ground cell, aiming at the center of the cell.

Once all the rays have been created, SamsaraLight performs the ray casting as described in Courbaud et al. (2003). For each ray, its energy is attenuated when it crosses a crown. The proportion of energy transmitted follows the formulation of the Beer-Lambert law:

$$\begin{aligned} I_T = I_0 e^{-k \times \Omega \times LAD \times l_p} \end{aligned}$$

(4)

where $l_p$ is the path length of the ray in the crown and $I_0$ is the energy of the beam before it intercepts the crown. Then, the energy absorbed by a crown $I_A$ is the complement of the transmitted energy:

$$\begin{aligned} I_A = I_0 - I_T \end{aligned}$$

(5)

Note that SamsaraLight does not take directly into account the reflection of rays—which causes a loss of energy in the sky and a reabsorption of the energy reflected on the ground. These phenomena are taken into account when calculating the attenuation coefficient (see Appendix Sect. 2).

After interception by a crown, the ray continues its course until it reaches either a new crown or a ground cell to which the remaining energy is transmitted. A proportion of absorbed radiation $\epsilon$ is uniformly removed from soil cells to represent the radiation interception from trunks, assuming a random arrangement of trees:

$$\begin{aligned} \epsilon = 1 - exp\left( - \frac{ \sum _i TS_i }{ S } \right) \end{aligned}$$

(6)

where S is the stand area and $\sum _i TS_i$ is the sum of the trunk shade surface of individual trees. $TS_i$ depends on the DBH and height of each tree i (supposing a cylindrical shape of the trunk), as well as on the hourly sun angle $\beta (h)$:

$$\begin{aligned} TS_i = DBH_i \times \frac{height_i}{tan(\beta (h))} \end{aligned}$$

(7)

At the end of the ray casting, the model integrates for each crown and soil cell the amount of direct and diffuse energy received over a year.

Computation of hourly absorbed energy

The hourly absorbed radiation of any element is then computed using the ray casting on the one hand and the hourly incident radiation on the other hand.

For each absorbing element i (a soil cell or a tree crown) and for each type of radiation (direct/diffuse, PAR/NIR), the energy it absorbs at hourly scale is given by the hourly incident radiation gr(h) and the fraction of energy absorbed annually by this element, $I_{Ay}(i)$, divided by the total energy absorbed by all elements j over the year:

$$\begin{aligned} I_A(h,i) = gr(h) \times \frac{I_{Ay}(i)}{ \sum _j{I_{Ay}(j)}} \end{aligned}$$

(8)

The value of $I_A(h,i)$ has then to be amended because the ray casting uses values of LAD that assume trees are at their period of full vegetation. A surplus of energy is then removed afterward from each tree according to their daily level of leaf development. This surplus is redistributed into other trees and soil cells, as described in Sect. 4.

Distribution among layers

Within a real-life tree, some leaves can receive a large amount of light—which leads to a saturation of the photosynthesis capacities—while others are in the shade. The saturation phenomenon (and more generally the concavity of the absorbed light-photosynthesis relationship) forbids calculating photosynthesis by considering an average level of light absorption for the whole canopy: this would bias upwards the estimation of photosynthesis (Leuning et al., 1995). In CASTANEA, the energy absorbed by the canopy is therefore distributed into five layers of leaves, in which the absorbed energy is assumed to be relatively homogeneous. The layers are themselves divided between leaves under direct light (called sun leaves) and leaves in the shade. The distribution of energy into the different layers is described in the “Distribution of radiations into canopy layers and into sun” section in Appendix 2.

2.2 Data set

To evaluate the simulations, we used an existing data set (GMAP forest plot design, Jourdan et al., 2019, 2020) composed of 39 beech, fir and beech-fir plots sampled between 2014 and 2016. Plots are distributed on three sites from the French pre-Alps (Bauges, Ventoux, Vercors), which are described in Table 2. They consist in a 10-m radius area in which the position, height, crown base height, age, diameter and species of each tree with a DBH greater than 7.5 cm were collected once.

Out of 1177 stems, 731 were cored to assess the growth dynamics over the 18-year period 1996–2013 (Jourdan et al., 2019). Growth of non-cored stems was inferred on the assumption that basal area increment over basal area was constant for a given species and site. To be comparable with the model output, basal area increments were converted into wood volume increments. To do that, we inferred past tree heights by using values of past DBH and the relationship between measured height and DBH. Past DBH were reconstructed using basal area increments and measured DBH. Then, a model was fitted on trees of the same species and site to evaluate the relationship between measured height and DBH (see Appendix 1). This model was used to compute past height based on reconstructed past DBH.

Wood volume increments were computed by multiplying each tree basal area increment with its inferred past height and $\Phi$, a form factor coefficients which takes into account the non-cylindrical shape of the trunks (Deleuze et al., 2014). On the one hand, PDG-Arena was evaluated using wood volume increments at individual scale. On the other hand, we used the wood volume increments at stand scale to evaluate both PDG-Arena and CASTANEA.

Hourly climate data (temperature, global radiation, wind speed, precipitation and relative humidity) were obtained from the 8-km scale SAFRAN reanalysis data set (Vidal et al., 2010) for the three sites and temperatures were adapted to each stand altitude using an adjustment of 0.6 °C per 100 m (Rolland, 2003). Soil texture, depth and stone content were obtained for every stand (data from one of the authors).

The LAI of the stands were retrieved using each plot coordinates and the 1-km resolution SPOT/PROBA-V remote sensing data set (Baret et al., 2013). We computed the average value of the yearly maximum LAI observed over the 1999–2013 period.

Table 2 Characteristics of the stands used to evaluate the model

2.3 Simulation plan

Using field inventories, we generated three sets of virtual inventories denoted RM, R, and O, following three levels of simplification compared to reality. The first set represents regularized monospecific inventories (RM): for each species of each stand, we generated a new inventory with equally spaced trees of the same species, age, diameter, and height. For mixed stands, the simulation results using RM inventories were assembled relatively to the proportion of each species basal area. RM inventories can then be used to simulate the growth of multispecific stands while ignoring species interactions. The second set represents regularized inventories (R), in which trees of different species can coexist but trees of the same species share the same age, diameter and height. Trees in R inventories are regularly spaced in a random order, independently of the species. Lastly, original inventories (O) include the information of the real life data set, that is species, position, diameter, and height of every individual trees. For each type of inventories representing the same stand (regularized or not, with or without species interactions), the mean quadratic diameter, volume per tree and tree age per species and the basal area were conserved.

CASTANEA was used as a reference model to evaluate the performance enhancement brought by PDG-Arena. We used RM inventories for CASTANEA’s stand-scale simulations. It is to be noted that, contrary to PDG-Arena, CASTANEA does not account for the stand slope. Therefore, when comparing CASTANEA and PDG-Arena results (Sect. 3.1), the slope was put to zero in PDG-Arena inventories. In the other situations (Sects. 3.2 and 3.3), the slopes of the inventories simulated using PDG-Arena were those of the field data.

To sum up, we simulated the growth of 39 stands over the 18-year period 1996–2013, considering four situations: RM, R, and O inventories with PDG-Arena and RM inventories with CASTANEA. Tree reproduction and intraspecific diversity, which are characteristics of PDG and therefore PDG-Arena, were switched off for these simulations.

2.4 Model evaluation

To evaluate the similarity between each modeling situation, we used the gross primary production (GPP) as CASTANEA and PDG-Arena are carbon-based models. We computed the coefficient of correlation (r, from − 1 to 1) for the simulated GPP per stand between the four situations (Sect. 2.3).

To evaluate the performance of the models against field measurements, we used the simulated wood volume increment per stand. We computed the Mean Absolute Percentage Error (MAPE) and the coefficient of determination (r², from 0 to 1) between simulations and measurements. A low MAPE indicates that simulated wood production is on average close to measured production. An r² close to 1 shows a good capacity of the model to predict stand production variability. Additionally, PDG-Arena with O inventories was evaluated at the individual tree scale, by computing the r² and MAPE of the simulated versus measured wood volume increment per tree for each group of the same site, type of stand (beech, fir of mixed) and species.

Lastly, we computed the net mixing effect (NME) to assess the extent of the simulated physiological processes that can solely be attributed to species mixing. Following the computation of the net biodiversity effect by Loreau (2010), we defined the NME as the difference for a variable between its observed value in mixed stands and its predicted value based on the hypothesis that there is no complementarity effect between species. Here, we compared the value of a simulated variable with PDG-Arena using the R and RM inventories (i.e., with and without species interactions). NME was evaluated on GPP, canopy absorbance (including PAR and NIR domains), transpiration rate and maximum water shortage (defined as the maximum difference reached during simulation between the current and full useful reserve, in mm). NME was tested against the null hypothesis using a two-sided Wilcoxon signed rank test.

3 Results

3.1 Comparison of PDG-Arena and CASTANEA

Using regular and monospecific inventories (RM), CASTANEA and PDG-Arena showed similar predictions for the stand-level GPP, with a coefficient of correlation at 99.8%. However, the GPP simulated by PDG-Arena was in average 4.2% larger than that by CASTANEA (Fig. 3). As shown in Table 3, which compares the 4 modeling situations based on the coefficient of correlation, simulations from PDG-Arena was closer to those of CASTANEA when using regularized inventories (R) on the one hand and when using regularized monospecific inventories (RM) on the other hand.

Fig. 3

Gross primary production (GPP) per stand simulated by PDG-Arena and CASTANEA. Regularized monospecific inventories (RM) were used. r is the correlation coefficient. Plots are represented as squares, circles and triangles for, respectively, Bauges, Vercors and Ventoux sites. They are colored in red, green and blue for, respectively, mixed, monospecific beech, and monospecific fir plots

Table 3 Matrix of similarity between simulated GPP from CASTANEA and PDG-Arena

3.2 Model performance

The simulated versus measured stand wood volume increment for the 39 stands are reported for the CASTANEA model using RM inventories and for the PDG-Arena model using O inventories in Appendix 3. Two fir stands from the Bauges site, denoted haut_sp_2 and bas_sp_4, stand out from the point cloud, with measured growths of 1995 and 1562 cm³m⁻²year⁻¹, while the simulated growth did not exceed 973 cm³m⁻²year⁻¹ for CASTANEA and PDG-Arena. In addition, simulations using values of LAI measured in 2022 using Terrestrial Laser Scanning (unpublished data from one of the authors) were performed and showed the same discrepancy with growth measurements for these two stands, indicating that the LAI source of measurement does not explain the poor performance on outliers. As the inclusion of these two stands in the analysis affects the overall results, they were discarded from the following analysis (see Appendix 3 Table 6 for the performance analysis that includes all stands).

Table 4 Evaluation of the performances of PDG-Arena and CASTANEA on the 37 stands

Simulation performances of CASTANEA and PDG-Arena against measured wood volume increments per stand are reported in Table 4. The MAPE was close between models and types of inventories, ranging from 30.1 to 33.1% in mixed stands, 53.9 to 57.9% in beech stands and 29.6 to 33.7% in fir stands. Considering the 37 stands, performances were close between CASTANEA and PDG-Arena on comparable inventories, i.e., RM inventories, with a slight advantage for PDG-Arena (r² 32.1% vs 29.5%). Using O inventories, PDG-Arena performed better than CASTANEA on RM inventories (r² 34.2 vs 29.5%).

Activation of species interactions in PDG-Arena (R vs RM inventories) slightly decreased the performance for mixed stands (r² 36.3% vs 37.6%, MAPE 33.1% vs 30.7%). Using original instead of regularized inventories (O vs R), PDG-Arena displayed an improved performance on mixed (r² 40.5 vs 36.3%, MAPE 31.5 vs 33.1%) and beech (r² 38.3 vs 24.7%, MAPE 53.9 vs 57.9%) stands but a lower performance on fir stands (r² 39.8 vs 50.1%, MAPE 39.8 vs 33.0%).

Appendix 3 Fig. 8 shows the simulated versus measured wood volume increment per year during the period 1996–2013 at tree scale using PDG-Arena and original inventories (O). r² ranged from 20 to 64% depending on the set of trees, with a mean at 47%. MAPE ranged from 50 to 146%, with a mean of 82% (Table 7).

3.3 Mixing and structure effects

GPP and canopy absorbance simulated by PDG-Arena in mixed stands are represented in Fig. 4 for RM, R and O inventories. Additionally, Appendix 3 Fig. 9 shows the yearly transpiration rate and maximum water shortage. Comparison of simulations with R and RM inventories showed a positive net mixing effect of 5.5% on GPP (1665 vs 1578 gC m⁻²year⁻¹; p-value < 0.001), of 11.1% on canopy absorbance (0.452 vs 0.407; p-value < 0.001), of 15.8% on canopy transpiration (234 vs 202 mm year⁻¹; p-value < 0.001) and of 13.7% on maximum water shortage (92.5 vs 81.3 mm; p-value < 0.001).

The structure effect (evaluated by comparing O and R inventories on all 39 stands, not shown here) decreased GPP by 3.7% (1603 vs 1665 gC m⁻².year⁻¹; p-value < 0.001) and the canopy absorbance by 5.2% (0.428 vs 0.452; p-value < 0.001). Transpiration showed a decrease of 3.2% (226 vs 234 mm; p-value < 0.001) and maximum water shortage a decrease of 1.9% (90.8 vs 92.51 mm; p-value < 0.05).

Fig. 4

Gross primary production (GPP) and canopy absorbance simulated by PDG-Arena for 13 mixed stands. Three types of inventories were used: regularized monospecific inventories (RM, in very light blue), regularized inventories with species interactions (R, in light blue) and original inventories (O, in dark blue). Two-sided Wilcoxon signed rank test was used (ns: not significant, *: p-value < 0.05, ***: p-value < 0.001)

4 Discussion

Given the paucity of forest growth models simulating ecophysiological processes at the individual tree scale, we developed the individual-based model PDG-Arena from the stand-scale model CASTANEA in order to simulate the carbon, water, and irradiance interception dynamics of mixed forests. PDG-Arena was built with the idea of observing and understanding the properties that emerge in multispecific stands, by integrating tree-level competition and without assuming the occurrence of positive interactions between heterospecific trees. It uses on the one hand a physiological model parameterized for monospecific stands and on the other hand an individual scale structure that allows trees to interact - the interaction being more or less competitive depending on the functional traits of the individuals and species.

We showed that PDG-Arena was able to reproduce the behavior of CASTANEA when simulating regularized inventories with no species interactions. Thus, the increase in complexity of PDG-Arena, required in order to simulate the functioning and interactions of distinct trees, was not at the cost of decreased performance at stand scale. Even when using original inventories (i.e., integrating the diversity in structure and species), the stand scale results of PDG-Arena were highly correlated to those of CASTANEA. This is explained by the fact that both models are based on LAI, which remains identical for each stand between simulations. Still, PDG-Arena, in comparison to CASTANEA, showed an improved performance when compared to measurements, in particular on beech (r² + 15.4 percentage points) and mixed stands (r² + 4.2 percentage points). As shown by the simulations using different types of inventories, the improvement in simulating stand growth is largely explained by the use of original stand structures, letting PDG-Arena simulate the growth of trees of various sizes.

At the individual scale, PDG-Arena explained half of the variability of tree growth, showing that it can capture the competitive status of each tree based on their leaf surface, height, and position. However, the mean absolute error was often large and systematic, indicating that the model would have to be calibrated more finely on each sites if it were to be used for predictions.

Interestingly, a positive and significant net mixing effect was observed in PDG-Arena simulations on gross primary productivity by comparing simulations with interacting species to equivalent simulations with species in isolation. The simulated overyielding can be attributed to an improvement of canopy absorbance due to species mixing (Fig. 4). LAI being equal between each inventory for the same stand, the increased radiation absorption is hence explained by a greater occupation of the aerial space due to species interactions. This effect, known as canopy packing, has been observed on a variety of mixed forests across Europe (Jucker et al., 2015; Pretzsch, 2019). Canopy packing is commonly decomposed into two processes: the phenotypic plasticity of the shape and size of crowns and the vertical stratification (i.e., the occupation by crowns of different vertical strata). Although it is likely to play a role in the functioning of mixed stands (Pretzsch, 2019; Dieler and Pretzsch, 2013), phenotypic plasticity is not yet implemented in PDG-Arena. Thus, our model can only simulate the vertical stratification of crowns, but not their morphological adaptation to their local competitor (see, for example, Jonard et al., 2020 and Morin et al., 2021), potentially leading to an underestimation of overyielding.

The observed overyielding in the French National Forest Inventory for beech-fir mixtures (20%, Toïgo et al., 2015) is larger than the one we simulated. In addition to canopy packing, the real-life overyielding in mixed stands can also be explained by reduced competition for nutrients. Indeed, nutrient content in above-ground biomass and the nitrogen concentration of leaves are likely to be increased by species mixing (Richards et al., 2010). However, competition for nutrients was not integrated in PDG-Arena since its main objective was to build an individual-based model upon the physiological processes that already exist in CASTANEA.

In addition, species mixing increased the yearly water shortage due to increased transpiration and irradiance interception at equivalent LAI (Appendix 3 Fig. 9). This confirms the idea that the diversity-functioning relationship in forests has differing effects depending on the resource considered (Forrester, 2014). According to our simulations, promoting diverse stands could maximize light interception and growth but would also increase transpiration, which would be detrimental in sites with limited water reserves. In reality, an increase in water use in mixed stands could be counter-balanced by a reduced competition for water between trees of different species (Schume et al., 2004). Although an interspecific differentiation between the water uptake depth has been observed for some species (Schwendenmann et al., 2015), our model cannot simulate this mechanism yet. A comprehensive knowledge of each species water uptake depth is still in construction but could be integrated in process-based models in the near future (Bachofen et al., 2024). Concerning the horizontal distance of tree water uptake, little data exist at the moment. The assumption of a horizontally homogeneous water uptake in our model is justified by the small surface area of the simulated plot.

One limit of this study was the nature of the data used to evaluate the model. Tree growth is an integrative measure that results from carbon assimilation, water uptake and irradiance interception, whereas CASTANEA was calibrated using CO₂ fluxes (Dufrêne et al., 2005). Moreover, the modeling of carbon allocation, which plays a decisive role in simulating wood growth, is a potential source of error (Davi et al., 2009; Merganičová et al., 2019). Additionally, climate was parameterized at the site scale using an 8-km resolution data set instead of at the stand scale, although climatic variables such as precipitation could vary between stands due to local topography.

5 Conclusion

The new individual-based model PDG-Arena we developed is able to simulate the interactions between trees in monospecific and mixed stands and predict their productivity based on an explicit tree inventory. Compared to CASTANEA, PDG-Arena showed improved predictive capability for beech and mixed beech-fir forests. The model can simulate the growth of small-sized stands (less than 1 ha), of regular or irregular structure, and composed of trees of similar or different species, given that the species ecophysiological properties are parameterized in CASTANEA. So far, CASTANEA has been parameterized for common beech (Dufrêne et al., 2005), silver fir (Davi and Cailleret, 2017), Douglas fir (Brèteau-Amores et al., 2019), Atlas cedar (Journé et al., 2021), Scotch pine and holm oak (Delpierre et al., 2012), maritime pine (Petit-Cailleux et al., 2021), and spruce, sessile oak and pedunculate oak (Guillemot et al., 2017). As PDG-Arena simulates the competition for water and light between trees with no preconceived ideas about the direction of interspecific interaction (from competition to complementarity), it can be used to test specific hypotheses about mixed forests and better understand the diversity-functioning relationship in forests under contrasted scenarios. For example, the model could be used to explore the following open questions, keeping in mind that the answers are largely species-specific and environment-dependent (Ratcliffe et al., 2015; Forrester et al., 2017): is overyielding more likely to occur in less productive sites (Toïgo et al., 2015)? Can overyielding increase water stress in mixed stands (Forrester et al., 2016)? Are mixed stands more resilient to drought (Grossiord, 2019)? Lastly, being built on the basis of a physio-demo-genetics model (Oddou-Muratorio and Davi, 2014), PDG-Arena could be used to evaluate the evolutionary dynamics of functional traits of a population under various biotic (stand composition, density, and structure) and abiotic (soil, climate) constraints, as intraspecific diversity is a major adaptive force in natural tree populations (Lefèvre et al., 2014; Oddou-Muratorio et al., 2020; Fady et al., 2020).

Appendices

Appendix 1 Height-diameter relationship

For each group of trees of the same species and site, a linear model (Eq. 9) was fitted on the logarithms of their measured height (in m) and DBH (in cm) as shown in Fig. 5. The slope and intercept parameter a and b as well as the coefficients of determination r² are shown in Table 5 for each group.

$$\begin{aligned} log(height) = a \times log10(DBH) + b \end{aligned}$$

(9)

Fig. 5

Relationship between measured height and DBH for every sampled trees. The red line indicates the linear model fitted on logarithmic values for each group of trees per site and plot composition. $r^2$ is the coefficient of determination

Table 5 Parameters of the height-DBH model described in Eq. 9

Appendix 2 Supplementary description of PDG-Arena

1.1 Computing leaf mass per area

The leaf mass per area (LMA) is a leaf-level trait defined as the mass per unit area of leaves (g m⁻²). LMA varies both in time (during leaf growth) and in space (depending on the leaf competitive status). There is therefore an exponentially decreasing distribution of LMA across the canopy from top to bottom. This section describes how the spatial variation of LMA is accounted in PDG-Arena.

In the CASTANEA model, which assumes that the stand is homogeneous and monospecific, the LMA profile follows a exponentially decreasing function (Davi et al., 2008):

$$\begin{aligned} LMA(LAI_{above}) = LMA_0 \times e^{-kLMA \times LAI_{above} } \end{aligned}$$

(10)

$LAI_{above}$ is the Leaf Area Index that accounts only for the leaves that are above the considered leaf. $LMA_0$ and kLMA depend on the species and describe the decrease in LMA within the canopy, which is related to the decrease in local irradiance within the canopy. Then, an average value for LMA within a layer is obtained by integrating $LMA(LAI_{above})$ within the layer’s vertical boundaries.

In the case of PDG-Arena, the canopy is more structurally complex than in CASTANEA and can include several species. In the new model, LMA at a given position in a tree is computed taking into account all the trees in the plot, and using the same formula as in Eq. 10. $LAI_{above}$ is computed by counting only the leaves of the canopy that are located above the considered leaf. It should be noted that the model is not completely accurate given that the parameter kLMA and $LMA_0$ are those of the species of the considered leaf, although the leaves taken into account in $LAI_{above}$ may be from a different species. However, this method does represent the phenomenon of light attenuation which is specific to each individual.

1.2 Estimation of the attenuation coefficient with reverse-engineering

In order to estimate the value of the attenuation coefficients of each species in PDG-Arena, a preliminary simulation is carried out following the CASTANEA model to take advantage of SAIL, its radiation sub-model (Dufrêne et al., 2005). The preliminary simulation is performed for each species on a monospecific and regularized inventory (RM inventory, see Sect. 2.3). We define the attenuation coefficient $k_1$ at a given time as a function of the incident energy $I_0$, the energy transmitted by the vegetation $I_t$, and the Leaf Area Index LAI, following a Beer-Lambert model:

$$\begin{aligned} I_t = I_0 exp^{-k_1 \times LAI } \end{aligned}$$

(11)

which is equivalent to:

$$\begin{aligned} k_1 = \frac{1}{LAI} \times log\left( \frac{I_0}{I_t} \right) \end{aligned}$$

(12)

where $I_t$ is defined at any time as the difference between the incident energy and the energy absorbed by the vegetation.

The coefficient of attenuation which is used in SamsaraLight, denoted $k_2$, is not of the same nature as $k_1$. Indeed, in Eq. 11, we multiply $k_1$ by the LAI (considering an infinite, horizontally homogeneous, leaf layer) while SamsaraLight multiplies $k_2$ to the Leaf Area Density LAD and the beam path length within a finite, volumetric crown (see Eq. 4). Then, to go from one to the other, we must multiply $k_1$ by $sin(\beta )$ (with $\beta$ the angle of height of the sun):

$$\begin{aligned} k_2 = sin(\beta ) \times k_1 = sin(\beta ) \times \frac{1}{LAI} \times log\left( \frac{I_0}{I_t} \right) \end{aligned}$$

(13)

Coefficient $k_2$ depends on the height of the sun, but also on the frequency domain of the radiation. Indeed, the attenuation coefficient takes into account both the extinction of the rays (defined by the leaf and crown geometry) and the absorption by the leaves which depends on the irradiance frequency. In the following calculations, we distinguish the PAR (photosynthetically active radiation) domain and the NIR (near infrared radiation) domain. It is assumed that these two domains represent the bulk of the incident radiation. To sum up, the attenuation coefficient should depend on the species (leaf angle distribution and absorbance rate), the type of radiation (PAR/NIR, direct/diffuse) and the sun height angle ($\beta$), which specific to each simulated ray.

Based on the results of the preliminary CASTANEA simulation, which executes a radiation balance using the SAIL model, we can infer the value of the attenuation coefficients of the plot for direct and diffuse radiations using Eq. 13. In the preliminary simulation, we know for direct rays the value of the height angle $\beta$ at any hour. For diffuse rays, by definition $\beta$ takes every value between 0 and $\pi / 2$ at any hour, so we can not use the height angle information.

Direct rays

For direct radiation, we estimate an attenuation coefficient for each species by discriminating the PAR and NIR and defining 20 classes of attenuation coefficients corresponding to classes of the height angle $\beta$, equally distributed between 0 and $\pi / 2$. For each class i of $\beta$, we performed an average on the attenuation coefficients observed during the preliminary simulation for direct radiations:

$$\begin{aligned} k_{dir}(i) = \sum \limits _{h_i} \left[ sin(\beta _i) \times \frac{1}{LAI(h_i)} \times log \left( \frac{I_{0dir}(h_i)}{I_{tdir}(h_i)} \right) \right] \times \frac{1}{n(h_i)} \end{aligned}$$

(14)

where $k_{dir}(i)$ is the mean attenuation coefficient computed from the preliminary simulation results, for direct radiation of the height angle class i (which includes $n(h_i)$ hours). For a given hour of the year $h_i$, $LAI(h_i)$ is the daily Leaf Area Index of the plot, $I_{0dir}(h_i)$, is the incident direct energy and $I_{tdir}(h_i)$ is the direct energy transmitted through the canopy.

Diffuse radiation

For diffuse radiation, we discriminate the attenuation coefficient according to the species and radiation domain only. The attenuation coefficient for diffuse radiation $k_{dif}$ is assumed to be constant for any sun height angle. To switch from one formulation of the Beer-Lambert law to the other (Eq. 13), a value of $\beta$ is nevertheless needed. Considering that the distribution of diffuse rays along the $\beta$ height angles is uniform, we simplify the equation by using an average computation. Then, we use $\overline{sin(\beta )}$, the average of $sin(\beta )$ for $\beta$ going from 0 to $\pi /2$ (which is about 0.637). For a species and a radiative domain, we compute an average on every day of year of the observed attenuation coefficient during the preliminary simulation:

$$\begin{aligned} k_{dif} = \sum \limits _j \left[ \overline{sin(\beta )} \times \frac{1}{LAI(j)} \times log\left( \frac{I_{0dif}(j)}{I_{tdif}(j)} \right) \right] \times \frac{1}{365} \end{aligned}$$

(15)

with, for day j, LAI(j) the Leaf Area Index, $I_{0dif}(j)$ the incident diffuse energy and $I_{tdif}(j)$ the diffuse energy transmitted through canopy.

1.3 Distribution of radiations into canopy layers and between sun and shade leaves

In CASTANEA, the energy absorbed by the canopy is distributed into five layers of leaves, which are themselves divided into leaves in direct light (called sun leaves) and leaves in the shade. We present here how PDG-Arena operates the distribution of the absorbed energy by individual crowns.

Proportion of sun leaves of a tree

The proportion of sun leaves of a crown, i.e., of its leaves subjected to direct radiation, is given by a formula borrowed from the HETEROFOR model (Jonard et al., 2020). Two factors define the shading received by the leaves of a tree: on the one hand, the external shading provided by competing trees, giving the proportion of sun leaves $pSun_{ext}$; on the other hand, the internal shading provided by the own leaves of a tree, giving the proportion of sun leaves $pSun_{int}$.

The shading provided by the competitors is given by the ratio of the direct incident energy above the tree $I_{d0}(aboveTree)$ to the potential direct incident energy $I_{d0}(potential)$, which is computed by SamsaraLight by ignoring neighbors trees:

$$\begin{aligned} pSun_{ext} = \frac{I_{d0}(aboveTree)}{I_{d0}(potential)} \end{aligned}$$

(16)

The second quotient to be evaluated is the proportion of the leaves of the tree shaded by its own leaves. The shading by the leaves of the tree itself follows the same relationship as the direct radiation within the tree, that is to say a Beer-Lambert law:

$$\begin{aligned} pSun_{int}(l) = exp^{-k_{dir} l} \end{aligned}$$

(17)

where $pSun_{int}(l)$ is the proportion of sun leaves remaining after the radiation passes through the crown, with l the cumulative LAI encountered by the passing beam and $k_{dir}$ the tree extinction coefficient for direct PAR. The proportion of sun leaves at the crown entrance is supposed to be 1, ignoring leaves shaded by neighboring trees.

We can compute $LAI_{sun-int}$, the amount of leaves that are not shaded by leaves of the same tree. To do this, we need to integrate $pSun_{int}(l)$ for l ranging from 0 to LAI, the Leaf Area Index of the tree:

$$\begin{aligned} LAI_{sun-int} & = \int _{0}^{LAI} pSun_{int}(l) dl \nonumber \\ & = \int _{0}^{LAI} e^{- k_{dir} l} dl \nonumber \\ & = \left[ \frac{e^{- k_{dir} l}}{-k_{dir}} \right] _{0}^{LAI} \nonumber \\ & = \frac{1 - e^{- k_{dir} LAI}}{k_{dir}} \end{aligned}$$

(18)

Thus, $pSun_{int} = LAI_{sun-int} / LAI$ represents the proportion of leaf remaining in the light when shaded by the tree’s own leaves.

Finally, the proportion of sun leaves of a tree is $pSun_{tree} = pSun_{ext} \times pSun_{int}$.

Distribution of radiations by layer

If SamsaraLight allows us to know the amount of energy absorbed per tree according to each domain (PAR/NIR) and type of energy (direct/diffused), noted $E_{tree}$, it does not allow us to distribute this amount between layers, differentiating leaves with high interception and leaves with low interception. To do so, we firstly divide the leaf surface of a tree into n equal-sized layers, and we assume that the radiative characteristics are homogeneous within a layer. We define a distribution function $f_i$, that determines $E_i$, the amount of energy that is absorbed by layer i:

$$\begin{aligned} E_i = E_{tree} \times \frac{f_i}{\sum _{n} f_i } \end{aligned}$$

(19)

We assume that the distribution $f_i$ is affected by the radiation interception from leaf surface that is located above the layer (whether it belongs to other trees or to the same tree). Then, we define a simple stand-scale model that describes the level of energy transmitted through the stand using the Beer-Lambert law. At any level of height located under a quantity of leaves $LAI_{above}$, the proportion of radiation transmitted through these leaves is:

$$\begin{aligned} p_{rad}(LAI_{above}) = e^{-k_{st} \times LAI_{above}} \end{aligned}$$

(20)

with $k_{st}$ the stand level attenuation coefficient. $LAI_{above}$ is calculated by counting the amount of leaves above the leaf layer under consideration, knowing the position and shape of each individual. A homogeneous distribution of leaf density within each individual crown is assumed. We do not consider the plot slope in this calculation, i.e., only tree height defines whether the leaves of one tree are higher than those of another.

Finally, to calculate $f_i$, the fraction of energy absorbed by any layer i of a crown, we compute the average value of $p_{rad}$ inside the layer by integrating it within its boundaries $LAI_{above}(i-1)$ and $LAI_{above}(i)$:

$$\begin{aligned} \begin{array}{l} f_i = \frac{\int _{LAI_{above}(i-1)}^{LAI_{above}(i)} e^{-k_{st} LAI_{above}} dLAI_{above}}{LAI_{above}(i) - LAI_{above}(i-1)} \\ \Longleftrightarrow \\ f_i = \frac{e^{-k_{st} LAI_{above}(i-1)} - e^{-k_{st} LAI_{above}(i)}}{k_{st}(LAI_{above}(i) - LAI_{above}(i-1))} \end{array} \end{aligned}$$

(21)

The proportion $f_i$ is computed for each type of radiation (direct/diffuse and PAR/NIR).

1.4 Reduction of absorbed radiations in SamsaraLight

In SamsaraLight standard mode, the foliage is assumed to be at its maximum during the whole process. Thus, the energy absorbed by the trees when their leaf area is in reality lower must be revised downwards, especially for deciduous trees, which lose all their leaves in autumn. For each individual, a ratio depending on its LAI is computed each day to represent the evolution of its absorption level from 0 to 1. The level of absorption is supposed to follow the dynamic of the Beer-Lambert law:

$$\begin{aligned} ratio_{LAI} = \frac{1 - e^{- k \times LAI}}{1 - e^{- k \times LAI_{max}}} \end{aligned}$$

(22)

For each radiation domain, k is the attenuation coefficient of a tree and $ratio_{LAI}$ is applied to its absorbed energy to take off the surplus. Nevertheless, the removed energy must be redistributed, because if it had not been intercepted, this energy would have been distributed among the other absorbing elements (crowns or soil cells). At this point, it is no longer possible to know to which element the energy should be distributed. Then, the extracted energy is redistributed to all absorbing elements, proportionally to their level of absorbed energy (after reduction according to LAI), which represents their relative interception capacity.

Appendix 3 Supplementary results

Figures 6 and 7 show the simulated versus measured wood volume increment per stand for the 39 stands (including the outliers) using, respectively, the CASTANEA model with RM inventories and the PDG-Arena model with O inventories. Table 6 shows the performance of the models at stand scale based on the r² and MAPE coefficients, computed without discarding the two silver fir outlier stands.

Fig. 6

Simulated versus measured wood volume increment (WVI) for the 39 stands using the PDG-Arena model and original inventories. Plots are represented as squares, circles and triangles for, respectively, Bauges, Vercors, and Ventoux sites. They are colors in red, green, and blue for, respectively, mixed, monospecific beech, and monospecific fir plots. Labeled points are the outlier plots

Fig. 7

Simulated versus measured wood volume increment (WVI) for the 39 stands using the CASTANEA model and RM inventories. Plots are represented as squares, circles and triangles for, respectively, Bauges, Vercors, and Ventoux sites. They are colors in red, green, and blue for, respectively, mixed, monospecific beech, and monospecific fir plots. Labeled points are the outlier plots

Table 6 Evaluation of the performances of PDG-Arena and CASTANEA without discarding outliers

Figure 8 shows the simulated versus measured wood volume increment per tree for the 37 stands using the PDG-Arena model with O inventories. Table 7 shows the individual-scale performances in terms of r² and MAPE .

Table 7 Performance at the individual scale of the PDG-Arena model using original inventories (O)

Fig. 8

Simulated versus measured wood volume increment (WVI) per year for every cored trees. The PDG-Arena model with original inventories (O) was used. Beech trees are represented in red and fir trees in blue. Gray lines indicate the 1:1 line

Figure 9 shows the maximum water shortage during an average year (i.e., the maximum difference reached during a year between the current and full useful reserve, in mm) and yearly transpiration simulated by PDG-Arena for 13 mixed stands using RM, R, and O inventories.

Fig. 9

Maximum water shortage and yearly transpiration simulated by PDG-Arena for 13 mixed stands. Maximum water shortage is defined as the yearly maximum difference reached between the current and full useful reserve. Three types of inventories were used: regularized monospecific inventories (RM, in very light blue), regularized inventories with species interactions (R, in light blue) and original inventories (O, in dark blue). Two-sided Wilcoxon signed rank test was used (*: p-value < 0.05, ***: p-value < 0.001)