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

Camille Rouet, Hendrik Davi, Arsène Druel, Bruno Fady, Xavier Morin

► To cite this version:

Camille Rouet, Hendrik Davi, Arsène Druel, Bruno Fady, Xavier Morin. PDG-Arena: an ecophysiological model for characterizing tree-tree interactions in heterogeneous and mixed stands. Annals of Forest Science, 2025, 82 (1), pp.8. 10.1186/s13595-025-01277-8. hal-04973308v3

HAL Id: hal-04973308

https://hal.science/hal-04973308v3

Submitted on 3 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

RESEARCH PAPER

Open Access

<figure><img src="image_4.png" /></figure>

Check for
updates

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

Camille Rouet¹,²* ID, Hendrik Davi¹ ID, Arsène Druel¹ ID, Bruno Fady¹ ID and Xavier Morin³ ID

**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.

**Keywords** Process-based modeling, Mixed forest, Competition, Overyielding, Drought, French Alps

Handling editor: Erwin Dreyer.

This paper has been recommended as a pre-print by the Peer Community in Forest & Wood Sciences. The full recommendation and the reviews are available here: https://forestwoodsci.peercommunityin.org/articles/rec?id=147. The final text of the pre-print is available here: https://www.biorxiv.org/content/10.1101/2024.02.09.579667v4. The version published here has been slightly edited during the handling by the editorial board of Annals of Forest Science, a PCI-friendly journal.

contact.camille.rouet@pm.me

Full list of author information is available at the end of the article

## 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

<figure><img src="image_5.png" /></figure>

© The Author(s) 2025. **Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

<footnote>*Correspondence:
Camille Rouet</footnote>
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; Toigo 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.

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^{-2}$. 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^2$ 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).
<figure><img src="image_6.png" /><figcaption>**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</figcaption></figure>

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
<table><caption>**Table 1** Temporal and spatial scales of physical and physiological processes in PDG-Arena</caption><thead><tr><th></th><th>Tree level</th><th>Stand level</th></tr></thead><tbody><tr><td rowspan="4">Hourly level</td><td>Photosynthesis</td><td>Ray casting</td></tr><tr><td>Respiration</td><td>Soil evaporation</td></tr><tr><td>Crown transpiration</td><td></td></tr><tr><td>Crown evaporation</td><td></td></tr><tr><td rowspan="3">Daily level</td><td>Water interception</td><td>Water balance</td></tr><tr><td>Leaf phenology</td><td></td></tr><tr><td>Carbon allocation</td><td></td></tr><tr><td>Yearly level</td><td>Tree growth</td><td></td></tr></tbody></table>

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 omitted.

Considering direct precipitation, runoff along trunks and transmitted rainfall on the one hand, and tree evaporation, 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*;

<figure><img src="image_7.png" /><figcaption>**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</figcaption></figure>
* 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 Ω 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):

$$ \text{crownRadius} = \beta_{\text{crown}} + \alpha_{\text{crown}} \times \text{DBH} \quad (1) $$

$\alpha_{crown}$ and $\beta_{crown}$ are species dependent parameters estimated on site at Mont Ventoux (unpublished data from one of the authors). Ω 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 LAi) is determined using the stand leaf area at full vegetation (LAstand, 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 LAth:

$$ LA_i = LA_{stand} \times \frac{LA_{th}(DBH_i, species_i)}{\sum_{j}^{n} LA_{th}(DBH_j, species_j)} \quad (2) $$

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

$$ LA_{th}(DBH_i, species_i) = \beta_0(species_i) \times DBH^{\beta_1(species_i)} \quad (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:

$$ I_T = I_0 e^{-k \times \Omega \times LAD \times l_p} \quad (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:

$$ I_A = I_0 - I_T \quad (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:

$$ \epsilon = 1 - \exp\left(-\frac{\sum_i TS_i}{S}\right) \quad (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)$:

$$ TS_i = DBH_i \times \frac{height_i}{\tan(\beta(h))} \quad (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*<sub>Ay</sub>(*i*), divided by the total energy absorbed by all elements *j* over the year:

$$I_A(h, i) = gr(h) \times \frac{I_{Ay}(i)}{\sum_j I_{Ay}(j)} \qquad (8)$$

The value of *I*<sub>A</sub>(*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

<table><caption>**Table 2** Characteristics of the stands used to evaluate the model</caption><thead><tr><th>Group</th><th>N</th><th>Altitude</th><th>Mean DBH</th><th>Density</th><th>Basal area</th><th>% Beech</th><th>Mean age</th><th>LAI</th></tr></thead><tbody><tr><td>Bauges</td><td>10</td><td>1100 ± 101</td><td>28.7 ± 6.7</td><td>1030 ± 685</td><td>72 ± 14</td><td>53 ± 43</td><td>89 ± 16</td><td>5.6 ± 0.2</td></tr><tr><td>Vercors</td><td>14</td><td>1250 ± 101</td><td>32.3 ± 8.6</td><td>657 ± 275</td><td>56 ± 14</td><td>53 ± 38</td><td>118 ± 40</td><td>5.6 ± 0.3</td></tr><tr><td>Ventoux</td><td>15</td><td>1250 ± 126</td><td>22.1 ± 6.3</td><td>1450 ± 623</td><td>57 ± 13</td><td>50 ± 40</td><td>105 ± 47</td><td>3.2 ± 0.3</td></tr><tr><td>Mixed</td><td>13</td><td>1200 ± 131</td><td>26.2 ± 7.3</td><td>1080 ± 465</td><td>64 ± 13</td><td>46 ± 10</td><td>101 ± 29</td><td>4.7 ± 0.5</td></tr><tr><td>Beech</td><td>14</td><td>1230 ± 118</td><td>26.7 ± 10.3</td><td>1200 ± 794</td><td>56 ± 14</td><td>97 ± 5</td><td>119 ± 35</td><td>4.7 ± 1.2</td></tr><tr><td>Fir</td><td>12</td><td>1190 ± 139</td><td>29.8 ± 7.4</td><td>867 ± 578</td><td>62 ± 18</td><td>5 ± 7</td><td>94 ± 50</td><td>4.7 ± 1.3</td></tr><tr><td>All</td><td>39</td><td>1210 ± 126</td><td>27.5 ± 8.4</td><td>850 ± 632</td><td>60 ± 15</td><td>51 ± 39</td><td>105 ± 39</td><td>2.9 ± 1.2</td></tr></tbody></table>

Mean value and standard deviation for each site (Bauges, Ventoux, Vercors) and composition (Mixed, Beech, Fir) are shown for variables: number of stands, altitude (in m), mean diameter at breast height per stand (in cm), density (in stem ha⁻¹), basal area (in m² ha⁻¹), proportion of beech basal area (in %), mean age per stand (in years), Leaf Area Index (no unit)
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.

## 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.
<figure><img src="image_8.png" /><figcaption>**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</figcaption></figure>

<table><caption>**Table 3** Matrix of similarity between simulated GPP from CASTANEA and PDG-Arena</caption><thead><tr><th></th><th>CASTANEA (RM)</th><th>PDG-Arena (RM)</th><th>PDG-Arena (R)</th><th>PDG-Arena (O)</th></tr></thead><tbody><tr><td>CASTANEA (RM)</td><td>100.0</td><td>-</td><td>-</td><td>-</td></tr><tr><td>PDG-Arena (RM)</td><td>99.8</td><td>100.0</td><td>-</td><td>-</td></tr><tr><td>PDG-Arena (R)</td><td>99.3</td><td>99.5</td><td>100.0</td><td>-</td></tr><tr><td>PDG-Arena (O)</td><td>97.7</td><td>98.5</td><td>99.0</td><td>100.0</td></tr></tbody></table>

Different types of inventories are considered: “RM” (regularized and monospecific, i.e., without species interactions), “R” (regularized, but with species interactions), and “O” (original inventories). Similarity is expressed using the correlation coefficient (in %) of the simulated gross primary production for the 39 stands over the 1996–2013 period

## 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^3 m^{-2} year^{-1}$, while the simulated growth did not exceed 973 $cm^3 m^{-2} year^{-1}$ 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).

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^2$ 32.1% vs 29.5%). Using O inventories, PDG-Arena performed better than CASTANEA on RM inventories ($r^2$ 34.2 vs 29.5%).

Activation of species interactions in PDG-Arena (R vs RM inventories) slightly decreased the performance for mixed stands ($r^2$ 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^2$ 40.5 vs 36.3%, MAPE 31.5 vs 33.1%) and beech ($r^2$ 38.3 vs 24.7%, MAPE 53.9 vs 57.9%) stands but a lower performance on fir stands ($r^2$ 39.8 vs 50.1%, MAPE 39.8 vs 33.0%).
<table><caption>**Table 4** Evaluation of the performances of PDG-Arena and CASTANEA on the 37 stands</caption><thead><tr><th>Set</th><th>Model</th><th>Inventories</th><th>r<sup>2</sup></th><th>MAPE</th></tr></thead><tbody><tr><td rowspan="4">All stands</td><td>CASTANEA</td><td>RM</td><td>29.5</td><td>40.6</td></tr><tr><td>PDG-Arena</td><td>RM</td><td>32.1</td><td>40.5</td></tr><tr><td>PDG-Arena</td><td>R</td><td>32.5</td><td>41.8</td></tr><tr><td>PDG-Arena</td><td>O</td><td>34.2</td><td>40.4</td></tr><tr><td rowspan="4">Mixed</td><td>CASTANEA</td><td>RM</td><td>36.3</td><td>30.1</td></tr><tr><td>PDG-Arena</td><td>RM</td><td>37.6</td><td>30.7</td></tr><tr><td>PDG-Arena</td><td>R</td><td>36.3</td><td>33.1</td></tr><tr><td>PDG-Arena</td><td>O</td><td>40.5</td><td>31.5</td></tr><tr><td rowspan="4">Beech pure</td><td>CASTANEA</td><td>RM</td><td>22.9</td><td>55.3</td></tr><tr><td>PDG-Arena</td><td>RM</td><td>25.0</td><td>57.4</td></tr><tr><td>PDG-Arena</td><td>R</td><td>24.7</td><td>57.9</td></tr><tr><td>PDG-Arena</td><td>O</td><td>38.3</td><td>53.9</td></tr><tr><td rowspan="4">Fir pure</td><td>CASTANEA</td><td>RM</td><td>42.0</td><td>33.7</td></tr><tr><td>PDG-Arena</td><td>RM</td><td>51.9</td><td>29.6</td></tr><tr><td>PDG-Arena</td><td>R</td><td>50.1</td><td>30.4</td></tr><tr><td>PDG-Arena</td><td>O</td><td>39.8</td><td>33.0</td></tr></tbody></table>

Coefficient of determination ($r^2$, in %) and Mean Absolute Percentage Error (MAPE, in %) were computed for the simulated versus measured yearly wood volume increment per stand over the period 1996–2013. Inventories are characterized as: “RM” (regularized and monospecific, i.e., without species interactions), “R” (regularized, but with species interactions), and “O” (original inventories)

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^2$ 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 $g \cdot m^{-2} \cdot year^{-1}$; *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 \cdot year^{-1}$; *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 $g \cdot m^{-2} \cdot year^{-1}$; *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).

### 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

<figure><img src="image_9.png" /><figcaption>**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)</figcaption></figure>
## 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).

## 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.

$$ \log(\text{height}) = a \times \log10(\text{DBH}) + b \quad (9) $$

<figure><img src="image_10.png" /><figcaption>**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² is the coefficient of determination</figcaption></figure>
<table><caption>**Table 5** Parameters of the height-DBH model described in Eq. 9</caption><thead><tr><th>Site</th><th>Species</th><th>a</th><th>b</th><th>r<sup>2</sup></th></tr></thead><tbody><tr><td>Bauges</td><td>Beech</td><td>0.69</td><td>0.33</td><td>0.78</td></tr><tr><td>Bauges</td><td>Fir</td><td>0.81</td><td>0.065</td><td>0.86</td></tr><tr><td>Ventoux</td><td>Beech</td><td>0.62</td><td>0.31</td><td>0.62</td></tr><tr><td>Ventoux</td><td>Fir</td><td>0.72</td><td>0.097</td><td>0.81</td></tr><tr><td>Vercors</td><td>Beech</td><td>0.78</td><td>0.13</td><td>0.87</td></tr><tr><td>Vercors</td><td>Fir</td><td>0.83</td><td>0.033</td><td>0.90</td></tr></tbody></table>

## Appendix 2 Supplementary description of PDG-Arena

### 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 (gm⁻²). 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 an exponentially decreasing function (Davi et al., 2008):

$$LMA(LAI_{above}) = LMA_{0} \times e^{-kLMA \times LAI_{above}} \quad (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.

### 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:

$$I_{t} = I_{0} \exp^{-k_{1} \times LAI} \quad (11)$$

which is equivalent to:

$$k_{1} = \frac{1}{LAI} \times \log\left(\frac{I_{0}}{I_{t}}\right) \quad (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(β) (with β the angle of height of the sun):

$$k_{2} = \sin(\beta) \times k_{1} = \sin(\beta) \times \frac{1}{LAI} \times \log\left(\frac{I_{0}}{I_{t}}\right) \quad (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 (β), which specific to each simulated ray.

Based on the results of the preliminary CASTANEA simulation, which executes a radiation balance using the
SAII 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:

$$k_{dir}(i) = \sum_{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)} \quad (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 $\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:

$$k_{dif} = \sum_{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} \quad (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.

## 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:

$$pSun_{ext} = \frac{I_{d0}(aboveTree)}{I_{d0}(potential)} \quad (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:

$$pSun_{int}(l) = \exp^{-k_{dir}l} \quad (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 \\
&= \int_{0}^{LAI} e^{-k_{dir}l}dl \\
&= \left[ \frac{e^{-k_{dir}l}}{-k_{dir}} \right]_{0}^{LAI} \\
&= \frac{1 - e^{-k_{dir}LAI}}{k_{dir}}
\end{aligned} 
\quad (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:

$$ E_i = E_{tree} \times \frac{f_i}{\sum_n f_i} \quad (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:

$$ p_{rad}(LAI_{above}) = e^{-k_{st} \times LAI_{above}} \quad (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}
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)} \\
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{aligned} 
\quad (21) $$

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

## 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:

$$ ratio_{LAI} = \frac{1 - e^{-k \times LAI}}{1 - e^{-k \times LAI_{max}}} \quad (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 CAS-TANEA 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.
<figure><img src="image_11.png" /><figcaption>**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</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**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</figcaption></figure>
<table><caption>**Table 6** Evaluation of the performances of PDG-Arena and CASTANEA without discarding outliers</caption><thead><tr><th>Set</th><th>Model</th><th>Inventories</th><th>r<sup>2</sup></th><th>MAPE</th></tr></thead><tbody><tr><td rowspan="4">All stands</td><td>CASTANEA</td><td>RM</td><td>25.7</td><td>41.7</td></tr><tr><td>PDG-Arena</td><td>RM</td><td>26.5</td><td>41.6</td></tr><tr><td>PDG-Arena</td><td>R</td><td>26.4</td><td>42.8</td></tr><tr><td>PDG-Arena</td><td>O</td><td>24.0</td><td>41.7</td></tr><tr><td rowspan="4">Mixed</td><td>CASTANEA</td><td>RM</td><td>36.3</td><td>30.1</td></tr><tr><td>PDG-Arena</td><td>RM</td><td>37.6</td><td>30.7</td></tr><tr><td>PDG-Arena</td><td>R</td><td>36.3</td><td>33.1</td></tr><tr><td>PDG-Arena</td><td>O</td><td>40.5</td><td>31.5</td></tr><tr><td rowspan="4">Beech pure</td><td>CASTANEA</td><td>RM</td><td>22.9</td><td>55.3</td></tr><tr><td>PDG-Arena</td><td>RM</td><td>25.0</td><td>57.4</td></tr><tr><td>PDG-Arena</td><td>R</td><td>24.7</td><td>57.9</td></tr><tr><td>PDG-Arena</td><td>O</td><td>38.3</td><td>53.9</td></tr><tr><td rowspan="4">Fir pure</td><td>CASTANEA</td><td>RM</td><td>18.0</td><td>38.4</td></tr><tr><td>PDG-Arena</td><td>RM</td><td>24.8</td><td>34.9</td></tr><tr><td>PDG-Arena</td><td>R</td><td>23.7</td><td>35.6</td></tr><tr><td>PDG-Arena</td><td>O</td><td>19.1</td><td>38.6</td></tr></tbody></table>

Coefficient of determination ($r^2$, in %) and Mean Absolute Percentage Error (MAPE, in %) were computed for the simulated versus measured yearly wood volume increment per stand over the period 1996–2013. Inventories are characterized as: “RM” (regularized and monospecific, i.e., without species interactions), “R” (regularized, but with species interactions), and “O” (original inventories)

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^2$ and MAPE.

<table><caption>**Table 7** Performance at the individual scale of the PDG-Arena model using original inventories (O)</caption><thead><tr><th>Site</th><th>Stand type</th><th>Species</th><th>r<sup>2</sup></th><th>MAPE</th></tr></thead><tbody><tr><td>Bauges</td><td>Mixed</td><td>Beech</td><td>36</td><td>70</td></tr><tr><td>Bauges</td><td>Mixed</td><td>Fir</td><td>62</td><td>68</td></tr><tr><td>Bauges</td><td>Pure beech</td><td>Beech</td><td>64</td><td>63</td></tr><tr><td>Bauges</td><td>Pure fir</td><td>Fir</td><td>20</td><td>73</td></tr><tr><td>Ventoux</td><td>Mixed</td><td>Beech</td><td>40</td><td>95</td></tr><tr><td>Ventoux</td><td>Mixed</td><td>Fir</td><td>59</td><td>50</td></tr><tr><td>Ventoux</td><td>Pure beech</td><td>Beech</td><td>40</td><td>69</td></tr><tr><td>Ventoux</td><td>Pure fir</td><td>Fir</td><td>43</td><td>95</td></tr><tr><td>Vercors</td><td>Mixed</td><td>Beech</td><td>51</td><td>146</td></tr><tr><td>Vercors</td><td>Mixed</td><td>Fir</td><td>49</td><td>68</td></tr><tr><td>Vercors</td><td>Pure beech</td><td>Beech</td><td>51</td><td>115</td></tr><tr><td>Vercors</td><td>Pure fir</td><td>Fir</td><td>48</td><td>67</td></tr></tbody></table>

$r^2$ and MAPE, expressed in %, were computed on set of trees of the same site, type of stand and species

<figure><img src="image_13.png" /><figcaption>**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</figcaption></figure>
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.

<figure><img src="image_14.png" /><figcaption>**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)</figcaption></figure>

**Credits**

Figures were designed using images from flicon.com.

**License**

For the purpose of Open Access, the authors have applied a CC BY-NC 4.0 public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.

**Authors' contributions**

Camille Rouet: conceptualization, methodology, software, visualization, writing—original draft. Hendrik Davi: conceptualization, supervision. Arsène Druel: methodology, writing—review. Bruno Fady: project administration, supervision, writing—review. Xavier Morin: methodology, data curation, supervision, writing—review.

**Funding**

This work was financed by ADEME, the French Agency for Ecological Transition, and ONF, the French National Forests Office. The observation design used in this study is part of the GMAP network, partly funded by the OSU OREME (https://oreme.org/observation/foret/gmap/).

**Data availability**

The PDG-Arena model is part of the Capsis framework, which aims to facilitate collaborative and shared software development for forest science (Dufour-Kowalski et al., 2012). PDG-Arena is an extension of the Physio-Demo-Genetics model, whose code is only accessible to members of the Capsis community for testing purposes. PDG-Arena can be run through Capsis either by script or via a graphical interface, using a formalized inventory (describing each tree of

the plot, the plot and soil characteristics and physiological, demographic and genetic options), as well as a CASTANEA-adapted climate file.

A repository is made accessible on the Zenodo platform (Rouet et al., 2024b), containing:

* The script used for the generation of inventories;

* The inventories used in the simulations;

* Two data sets describing the LAI and soil textures used in the simulations;

* The results data set (including formatted data of simulations and measurements);

* The script used for the analysis of the result data set; Raw simulation results are accessible on demand.

**Declarations**

**Ethics approval and consent to participate**

Not applicable.

**Consent for publication**

Not applicable.

**Competing interests**

The authors of this publication declare that they have no conflicts of interest.

**Author details**

¹INRAE URFM, Domaine Saint Paul - Site Agroparc, Avignon CEDEX 9 84914, France. ²ADEME, 20 avenue du Grésillé, Angers CEDEX 1 49004, France. ³CEFE, CNRS, Univ. Montpellier, EPHE, IRD, 1919, route de Mende, Montpellier CEDEX 5 34293, France.
Received: 9 October 2024 Accepted: 14 January 2025
Published online: 25 February 2025

## References

Ammer C (2019) Diversity and forest productivity in a changing climate. New Phytol 221(1):50–66. https://doi.org/10.1111/nph.15263

Bachofen C, Tumber-Dávila SJ, Mackay DS, McDowell NG, Carminati A, Klein T, Stocker BD, Mencuccini M, Grossiord C (2024) Tree water uptake patterns across the globe. New Phytol 242(5):1891–1910. https://doi.org/10.1111/nph.19762

Barbosa LO, dos Santos JA, Gonçalves AFA, Campoe OC, Scolforo JRS, Scolforo HF (2023) Competition in forest plantations: Empirical and process-based modelling in pine and eucalypt plantations. Ecol Model 483:110410. https://doi.org/10.1016/j.ecolmodel.2023.110410

Baret F, Weiss M, Lacaze R, Camacho F, Makhmara H, Pacholcyzk P, Smets B (2013) GEOV1: LAI and FAPAR essential climate variables and FCOVER global time series capitalizing over existing products. Part 1: Principles of development and production. Remote Sens Environ 137:299–309. https://doi.org/10.1016/j.rse.2012.12.027

Bauhus J, Forrester DI, Pretzsch H (2017) From Observations to Evidence About Effects of Mixed-Species Stands. In: Pretzsch H, Forrester DI, Bauhus J (eds) Mixed-Species Forests: Ecology and Management, Springer, Berlin, Heidelberg, pp 27–71. https://doi.org/10.1007/978-3-662-54553-9_2

Bonan GB (2008) Forests and Climate Change: Forcings, Feedbacks, and the Climate Benefits of Forests. Science 320(5882):1444–1449. https://doi.org/10.1126/science.1155121

Brêteau-Amores S, Brunette M, Davi H (2019) An Economic Comparison of Adaptation Strategies Towards a Drought-induced Risk of Forest Decline. Ecol Econ 164:106294. https://doi.org/10.1016/j.ecolecon.2019.04.006

Cordonnier T, Smadi C, Kunstler G, Courbaud B (2019) Asymmetric competition, ontogenetic growth and size inequality drive the difference in productivity between two-strata and one-stratum forest stands. Theor Popul Biol 130:83–93. https://doi.org/10.1016/j.tpb.2019.07.001

Courbaud B, de Coligny F, Cordonnier T (2003) Simulating radiation distribution in a heterogeneous Norway spruce forest on a slope. Agric For Meteorol 116(1):1–18. https://doi.org/10.1016/S0168-1923(02)00254-X

Cuddington K, Fortin MJ, Gerber LR, Hastings A, Liebhold A, O'Connor M, Ray C (2013) Process-based models are required to manage ecological systems in a changing world. Ecosphere 4(2):art20. https://doi.org/10.1890/ES12-00178.1

Dănescu A, Albrecht AT, Bauhus J (2016) Structural diversity promotes productivity of mixed, uneven-aged forests in southwestern Germany. Oecologia 182(2):319–333. https://doi.org/10.1007/s00442-016-3623-4

Davi H, Cailleret M (2017) Assessing drought-driven mortality trees with physiological process-based models. Agric For Meteorol 232:279–290. https://doi.org/10.1016/j.agrformet.2016.08.019

Davi H, Baret F, Huc R, Dufrêne E (2008) Effect of thinning on LAI variance in heterogeneous forests. For Ecol Manag 256(5):890–899. https://doi.org/10.1016/j.foreco.2008.05.047

Davi H, Barbaroux C, Dufrêne E, François C, Montpied P, Bréda N, Badeck F (2008) Modelling leaf mass per area in forest canopy as affected by prevailing radiation conditions. Ecol Model 211(3):339–349. https://doi.org/10.1016/j.ecolmodel.2007.09.012

Davi H, Barbaroux C, Francois C, Dufrêne E (2009) The fundamental role of reserves and hydraulic constraints in predicting LAI and carbon allocation in forests. Agric For Meteorol 149(2):349–361. https://doi.org/10.1016/j.agrformet.2008.08.014

del Río M, Pretzsch H, Ruiz-Peinado R, Jactel H, Coll L, Löf M, Aldea J, Ammer C, Avdagić A, Barbeito I, Bielak K, Bravo F, Brazaitis G, Cerný J, Collet C, Condés S, Drössler L, Fabrika M, Heym M, Holm SO, Hylen G, Jansons A, Kurylyak V, Lombardi F, Matović B, Metslaid M, Motta R, Nord-Larsen T, Nothdurft A, den Ouden J, Pach M, Pardos M, Poeydebat C, Ponette Q, Pétor T, Reventlow DOJ, Sitko R, Sramek V, Steckel M, Svoboda M, Verheyen K, Vospernik S, Wolff B, Zlatanov T, Bravo-Oviedo A (2022) Emerging stability of forest productivity by mixing two species buffers temperature destabilizing effect. J App Ecol 59(11):2730–2741. https://doi.org/10.1111/1365-2664.14267

Deleuze C, Morneau F, Renaud J, Vivien Y, Rivoire M, Santenoise P, Longue-taud F, Mothe F, Hervé J, Vallet P (2014) Estimer le volume total d'un arbre, quelles que soient l'essence, la taille, la sylviculture, la station.

Rendez-vous Techniques de l'ONF 44:22-32. https://hal.science/hal-01143797v1

Delpierre N, Soudani K, François C, Le Maire G, Bernhofer C, Kutsch W, Misson L, Rambal S, Vesala T, Dufrêne E (2012) Quantifying the influence of climate and biological drivers on the interannual variability of carbon exchanges in European forests through process-based modelling. Agric For Meteorol 154–155:99–112. https://doi.org/10.1016/j.agrformet.2011.10.010

Dieler J, Pretzsch H (2013) Morphological plasticity of European beech (Fagus sylvatica L.) in pure and mixed-species stands. For Ecol Manag 295:97–108. https://doi.org/10.1016/j.foreco.2012.12.049

Dufour-Kowalski S, Courbaud B, Dreyfus P et al (2012) Capsis: an open software framework and community for forest growth modelling. Ann For Sci 69(2):221–233. https://doi.org/10.1007/s13595-011-0140-9

Dufrêne E, Davi H, François C, Maire GI, Dantec VL, Granier A (2005) Modelling carbon and water cycles in a beech forest. Part I: Model description and uncertainty analysis on modelled NEE. Ecol Model 185(2-4):407–436. https://doi.org/10.1016/j.ecolmodel.2005.01.00

Fady B, Aravanopoulos F, Benavides R, González-Martínez S, Grivet D, Lascoux M, Lindner M, Rellstab C, Valladares F, Vinceti B (2020) Genetics to the rescue: managing forests sustainably in a changing world. Tree Genet Genomes 16(6):80. https://doi.org/10.1007/s11295-020-01474-8

Fatichi S, Leuzinger S, Körner C (2014) Moving beyond photosynthesis: from carbon source to sink-driven vegetation modeling. New Phytol 201(4):1086–1095. https://doi.org/10.1111/nph.12614

Fontes L, Bontemps JD, Bugmann H, Mv Oijen, Gracia C, Kramer K, Lindner M, Rötzer T, Skovsgaard JP (2010) Models for supporting forest management in a changing environment. For Syst 19:8–29

Forrester DI (2014) The spatial and temporal dynamics of species interactions in mixed-species forests: From pattern to process. For Ecol Manag 312:282–292. https://doi.org/10.1016/j.foreco.2013.10.003

Forrester DI, Bauhus J (2016) A Review of Processes Behind Diversity-Productivity Relationships in Forests. Curr For Rep 2(1):45–61. https://doi.org/10.1007/s40725-016-0031-2

Forrester DI, Bonal D, Dawud S, Gessler A, Granier A, Pollastrini M, Grossiord C (2016) Drought responses by individual tree species are not often correlated with tree species diversity in European forests. J Appl Ecol 53(6):1725–1734. https://doi.org/10.1111/1365-2664.12745

Forrester DI, Tachauer IHH, Annighöfer P, Barbeito I, Pretzsch H, Ruiz-Peinado R, Stark H, Vacchiano G, Zlatanov T, Chakraborty T, Saha S, Sileshi GW (2017) Generalized biomass and leaf area allometric equations for European tree species incorporating stand structure, tree age and climate. For Ecol Manag 396:160–175. https://doi.org/10.1016/j.foreco.2017.04.011

Forrester DI, Ammer C, Annighöfer PJ, Avdagić A, Barbeito I, Bielak K, Brazaitis G, Coll L, del Río M, Drössler L, Heym M, Hurt V, Löf M, Matović B, Meloni F, den Ouden J, Pach M, Pereira MG, Ponette Q, Pretzsch H, Skrzyszewski J, Stojanović D, Svoboda M, Ruiz-Peinado R, Vacchiano G, Verheyen K, Zlatanov T, Bravo-Oviedo A (2017) Predicting the spatial and temporal dynamics of species interactions in Fagus sylvatica and Pinus sylvestris forests across Europe. For Ecol Manag 405:112–133. https://doi.org/10.1016/j.foreco.2017.09.029

Gonçalves AFA, Santos JAd, França LCdJ, Campoe OC, Altoé TF, Scolforo JRS (2021) Use of the process-based models in forest research: a bibliometric review. CERNE 27:e. https://doi.org/10.1590/0104776020217012769

González de Andrés E (2019) Interactions between Climate and Nutrient Cycles on Forest Response to Global Change: The Role of Mixed Forests. Forests 10(8):609. https://doi.org/10.3390/f10080609

Granier A, Bréda N, Biron P, Villette S (1999) A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands. Ecol Model 116(2):269–283. https://doi.org/10.1016/S0304-3800(98)00205-1

Grossiord C (2019) Having the right neighbors: how tree species diversity modulates drought impacts on forests. New Phytol 228(1):42–49. https://doi.org/10.1111/nph.15667

Grossiord C, Granier A, Ratcliffe S, Bouriaud O, Bruelheide H, Chečko E, Forrester DI, Dawud SM, Finér L, Pollastrini M, Scherer-Lorenzen M, Valladares F, Bonal D, Gessler A (2014) Tree diversity does not always improve resistance of forest ecosystems to drought. Proc Natl Acad Sci 111(41):14812–14815. https://doi.org/10.1073/pnas.1411970111

Guillemot J, Francois C, Hmimina G, Dufrêne E, Martin-StPaul NK, Soudani K, Marie G, Ourcival JM, Delpierre N (2017) Environmental control of carbon
allocation matters for modelling forest growth. New Phytol. 214:180-193.
https://doi.org/10.1111/nph.14320

Guillemot J, Kunz M, Schnabel F, Fichtner A, Madsen CP, Gebauer T, Härdtle W, von Oheimb G, Potvin C (2020) Neighbourhood-mediated shifts in tree biomass allocation drive overyielding in tropical species mixtures. New Phytol 228(4):1256-1268. https://doi.org/10.1111/nph.16722

Jonard M, André F, de Coligny F, de Wergifosse L, Beudez N, Davi H, Ligot G, Ponette Q, Vincke C (2020) HETEROFOR 1.0: a spatially explicit model for exploring the response of structurally complex forests to uncertain future conditions – Part 1: Carbon fluxes and tree dimensional growth. Geosci Model Dev 13(3):905-935. https://doi.org/10.5194/gmd-13-905-2020

Jourdan M, Lebourgeois F, Morin X (2019) The effect of tree diversity on the resistance and recovery of forest stands in the French Alps may depend on species differences in hydraulic features. For Ecol Manag 450:117486.
https://doi.org/10.1016/j.foreco.2019.117486

Jourdan M, Kunstler G, Morin X (2020) How neighbourhood interactions control the temporal stability and resilience to drought of trees in mountain forests. J Ecol 108(2):666-677. https://doi.org/10.1111/1365-2745.13294

Jourdan M, Cordonnier T, Dreyfus P, Riond C, de Coligny F, Morin X (2021) Managing mixed stands can mitigate severe climate change impacts on French alpine forests. Reg Environ Chang 21(3):78. https://doi.org/10.1007/s10113-021-01805-y

Journé V, Papaix J, Walker E, Courbet F, Lefèvre F, Oddou-Muratorio S, Davi H (2021). A hierarchical Bayesian model to investigate trade-offs between growth and reproduction in a long-lived plant. *BioRxiv pre-print*, https://doi.org/10.1101/2021.01.26.428205

Jucker T, Bouriaud O, Avacaritei D, Dănilă I, Duduman G, Valladares F, Coomes DA (2014) Competition for light and water play contrasting roles in driving diversity-productivity relationships in Iberian forests. J Ecol 102(5):1202-1213. https://doi.org/10.1111/1365-2745.12276

Jucker T, Bouriaud O, Coomes DA (2015) Crown plasticity enables trees to optimize canopy packing in mixed-species forests. Funct Ecol 29(8):1078-1086. https://doi.org/10.1111/1365-2435.12428

Korzukhin MD, Ter-Mikaelian MT, Wagner RG (1996) Process versus empirical models: which approach for forest ecosystem management? Can J For Res 26(5):879-887. https://doi.org/10.1139/x26-096

Lefèvre F, Boivin T, Bontemps A, Courbet F, Davi H, Durand-Gillmann M, Fady B, Gauzere J, Gidoin C, Karam MJ, Lalagüe H, Oddou-Muratorio S, Pichot C (2014) Considering evolutionary processes in adaptive forestry. Ann For Sci 71(7):723-739. https://doi.org/10.1007/s13595-013-0272-1

Leuning R, Kellner FM, Pury DGGD, Schulze ED (1995) Leaf nitrogen, photosynthesis, conductance and transpiration: scaling from leaves to canopies. Plant Cell Environ 18(10):1183-1200. https://doi.org/10.1111/j.1365-3040.1995.tb00628.x

Liang J, Crowther TW, Picard N, Wiser S, Zhou M, Alberti G, Schulze ED, McGuire AD, Bozzato F, Pretzsch H, de Miguel S, Paquette A, Hérault B, Scherer-Lorenzen M, Barrett CB, Glick HB, Hengeveld GM, Nabuurs GJ, Pfautsch S, Viana H, Vibrans AC, Ammer C, Schall P, Verbyla D, Tchebakova N, Fischer M, Watson JV, Chen HYH, Lei X, Schelhaas MJ, Lu H, Gianelle D, Parfenova EI, Salas C, Lee E, Lee B, Kim HS, Bruelheide H, Coomes DA, Piotto D, Sunderland T, Schmid B, Gourlet-Fleury S, Sonké B, Tavani R, Zhu J, Brandl S, Vayreda J, Kitahara F, Searle EB, Neldner VJ, Ngugi MR, Baraloto C, Frizzera L, Balazy R, Oleksyn J, Zawila-Niedźwiecki T, Bouriaud O, Bussotti F, Finér L, Jaroszewicz B, Jucker T, Valladares F, Jagodzinski AM, Peri PL, Gonmadje C, Marthy W, O'Brien T, Martin EH, Marshall AR, Rovero F, Bitariho R, Niklaus PA, Alvarez-Loayza P, Chamuya N, Valencia R, Mortier F, Wortel V, Engone-Obiang NL, Ferreira LV, Odeke DE, Vasquez RM, Lewis SL, Reich PB (2016) Positive biodiversity-productivity relationship predominant in global forests. Science 354(6309). https://doi.org/10.1126/science.aaf8957

Lindner M, Maroschek M, Netherer S et al (2010) Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. For Ecol Manag 259(4):698-709. https://doi.org/10.1016/j.foreco.2009.09.023

Loreau M (2010) CHAPTER 3. Biodiversity and Ecosystem Functioning. In: From Populations to Ecosystems, Princeton University Press, pp 56–78. https://doi.org/10.1515/9781400834167.56

Mas E, Cochard H, Deluigi J, Didion-Gency M, Martin-StPaul N, Morcillo L, Valadares F, Vilagrosa A, Grossiord C (2024) Interactions between beech and oak seedlings can modify the effects of hotter droughts and the onset of hydraulic failure. New Phytol 241(3):1021-1034. https://doi.org/10.1111/nph.19358

Merganičová K, Merganič J, Lehtonen A, Vacchiano G, Sever MZO, Augustynczik ALD, Grote R, Kyselová I, Mäkelä A, Yousefpour R, Krejza J, Collatti A, Reyer CPO (2019) Forest carbon allocation modelling under climate change. Tree Physiol 39(12):1937–1960. https://doi.org/10.1093/treephys/tpz105

Messier C, Bauhus J, Sousa-Silva R, Auge H, Baeten L, Barsoum N, Bruelheide H, Caldwell B, Cavender-Bares J, Dhiedt E, Eisenhauer N, Ganade G, Gravel D, Guillemot J, Hall JS, Hector A, Hérault B, Jactel H, Koricheva J, Kreft H, Mereu S, Muys B, Nock CA, Paquette A, Parker JD, Perring MP, Ponette Q, Potvin C, Reich PB, Scherer-Lorenzen M, Schnabel F, Verheyen K, Weih M, Wollni M, Zemp DC (2022) For the sake of resilience and multifunctionality, let's diversify planted forests! Conserv Lett 15(1):e12829. https://doi.org/10.1111/conl.12829

Metz J, Annighöfer P, Schall P, Zimmermann J, Kahl T, Schulze ED, Ammer C (2016) Site-adapted admixed tree species reduce drought susceptibility of mature European beech. Global Change Biol 22(2):903-920. https://doi.org/10.1111/gcb.13113

Monteith JL (1965) Evaporation and environment. Symp Soc Exp Biol 19:205-234.

Morin X, Fahse L, Scherer-Lorenzen M, Bugmann H (2011) Tree species richness promotes productivity in temperate forests through strong complementarity between species: Species richness promotes forest productivity. Ecol Lett 14(12):1211–1219. https://doi.org/10.1111/j.1461-0248.2011.01691.x

Morin X, Fahse L, de Mazancourt C, Scherer-Lorenzen M, Bugmann H (2014) Temporal stability in forest productivity increases with tree diversity due to asynchrony in species dynamics. Ecol Lett 17(12):1526–1535. https://doi.org/10.1111/ele.12357

Morin X, Bugmann H, Coligny FD, Martin-StPaul N, Cailleret M, Limousin JM, Ourcival JM, Prevosto B, Simioni G, Toigo M, Vennetier M, Catteau E, Guillemot J (2021) Beyond forest succession: A gap model to study ecosystem functioning and tree community composition under climate change. Funct Ecol. https://doi.org/10.1111/1365-2435.13760

Oddou-Muratorio S, Davi H (2014) Simulating local adaptation to climate of forest trees with a Physio-Demo-Genetics model. Evol Appl 7(4):453–467. https://doi.org/10.1111/eva.12143

Oddou-Muratorio S, Davi H, Lefèvre F (2020) Integrating evolutionary, demographic and ecophysiological processes to predict the adaptive dynamics of forest tree populations under global change. Tree Genet Genomes 16(5):67. https://doi.org/10.1007/s11295-020-01451-1

Pardos M, del Rio M, Pretzsch H, Jactel H, Bielak K, Bravo F, Brazaitis G, Defossez E, Engel M, Godvod K, Jacobs K, Jansone L, Jansons A, Morin X, Nothdurft A, Oreti L, Ponette Q, Pach M, Riofrío J, Ruiz-Peinado R, Tomao A, Uhl E, Calama R (2021) The greater resilience of mixed forests to drought mainly depends on their composition: Analysis along a climate gradient across Europe. For Ecol Manag 481:118687. https://doi.org/10.1016/j.foreco.2020.118687

Petit-Cailleux C, Davi H, Lefèvre F, Verkerk PJ, Fady B, Lindner M, Oddou-Muratorio S (2021) Tree mortality risks under climate change in Europe: assessment of silviculture practices and genetic conservation networks. Front Ecol Evol 9:582. https://doi.org/10.3389/fevo.2021.706414

Piotto D (2008) A meta-analysis comparing tree growth in monocultures and mixed plantations. For Ecol Manag 255(3):781-786. https://doi.org/10.1016/j.foreco.2007.09.065

Porté A, Bartelink HH (2002) Modelling mixed forest growth: a review of models for forest management. Ecol Model 150(1):141–188. https://doi.org/10.1016/S0304-3800(01)00476-8

Pretzsch H (2019) The Effect of Tree Crown Allometry on Community Dynamics in Mixed-Species Stands versus Monocultures. A Review and Perspectives for Modeling and Silvicultural Regulation. Forests 10(9):810. https://doi.org/10.3390/f10090810

Pretzsch H, Schütze G, Uhl E (2013) Resistance of European tree species to drought stress in mixed versus pure forests: evidence of stress release by inter-specific facilitation. Plant Biol 15(3):483–495. https://doi.org/10.1111/j.1438-8677.2012.00670.x

Pretzsch H, Forrester DI, Rötzer T (2015) Representation of species mixing in forest growth models. A review and perspective. Ecol Model 313:276–292. https://doi.org/10.1016/j.ecolmodel.2015.06.044

Pretzsch H, Forrester DI, Bauhus J (eds) (2017) Mixed-species forests: ecology and management. Springer, Berlin
Ratcliffe S, Holzwarth F, Nadrowski K, Levick S, Wirth C (2015) Tree neighbourhood matters – Tree species composition drives diversity–productivity patterns in a near-natural beech forest. For Ecol Manag 335:225–234. https://doi.org/10.1016/j.foreco.2014.09.032

Ratcliffe S, Liebergesell M, Ruiz-Benito P, González JM, Castañeda JMM, Kändler G, Lehtonen A, Dahlgren J, Kattge J, Peñuelas J, Zavala MA, Wirth C (2016) Modes of functional biodiversity control on tree productivity across the European continent. Glob Ecol Biogeogr 25(3):251–262. https://doi.org/10.1111/geb.12406

Reyer C (2015) Forest Productivity Under Environmental Change–a Review of Stand-Scale Modeling Studies. Curr For Rep 1(2):53–68. https://doi.org/10.1007/s40725-015-0009-5

Richards AE, Forrester DI, Bauhus J, Scherer-Lorenzen M (2010) The influence of mixed tree plantations on the nutrition of individual species: a review. Tree Physiol 30(9):1192–1208. https://doi.org/10.1093/treephys/tpq035

Rolland C (2003) Spatial and Seasonal Variations of Air Temperature Lapse Rates in Alpine Regions. J Clim 16(7):1032–1046. https://doi.org/10.1175/1520-0442(2003)016<1032:SASVOA>2.0.CO;2

Rouet C, Morin X, Druel A (2024b) Inputs, results data and analysis script for the evaluation of the PDG-Arena forest growth model on beech-fir stands. Data set: https://doi.org/10.5281/zenodo.12191049

Schume H, Jost G, Hager H (2004) Soil water depletion and recharge patterns in mixed and pure forest stands of European beech and Norway spruce. J Hydrol 289(1):258–274. https://doi.org/10.1016/j.jhydrol.2003.11.036

Schwendenmann L, Pendall E, Sanchez-Bragado R, Kunert N, Hölscher D (2015) Tree water uptake in a tropical plantation varying in tree diversity: interspecific differences, seasonal shifts and complementarity. Ecohydrology 8(1):1–12. https://doi.org/10.1002/eco.1479

Seynave I, Bailly A, Balandier P, Bontemps JD, Cailly P, Cordonnier T, Deleuze C, Dhôte JF, Ginisty C, Lebourgeois F, Merzeau D, Paillassa E, Perret S, Richter C, Meredieu C (2018) GIS Coop: networks of silvicultural trials for supporting forest management under changing environment. Ann For Sci 75(2):1–20. https://doi.org/10.1007/s13595-018-0692-z

Toïgo M, Vallet P, Perot T, Bontemps JD, Piedallu C, Courbaud B (2015) Overyielding in mixed forests decreases with site productivity. J Ecol 103(2):502–512. https://doi.org/10.1111/1365-2745.12353

Trogisch S, Liu X, Rutten G, Xue K, Bauhus J, Brose U, Bu W, Cesarz S, Chesters D, Connolly J, Cui X, Eisenhauer N, Guo L, Haider S, Härdtle W, Kunz M, Liu L, Ma Z, Neumann S, Sang W, Schuldt A, Tang Z, van Dam NM, von Oheimb G, Wang MQ, Wang S, Weinhold A, Wirth C, Wubet T, Xu X, Yang B, Zhang N, Zhu CD, Ma K, Wang Y, Bruelheide H (2021) The significance of tree-tree interactions for forest ecosystem functioning. Basic Appl Ecol 55:33–52. https://doi.org/10.1016/j.baae.2021.02.003

Trumbore S, Brando P, Hartmann H (2015) Forest health and global change. Science 349(6250):814–818. https://doi.org/10.1126/science.aac6759

van der Plas F, Manning P, Allan E, Scherer-Lorenzen M, Verheyen K, Wirth C, Zavala MA, Hector A, Ampoorter E, Baeten L, Barbaro L, Bauhus J, Benavides R, Benneter A, Berthold F, Bonal D, Bouriaud O, Bruelheide H, Bussotti F, Carnol M, Castagneyrol B, Charbonnier Y, Coomes D, Coppi A, Bastias CC, Muhie Dawud S, De Wandeler H, Domisch T, Finér L, Gessler A, Granier A, Grossiord C, Guyot V, Hättenschwiler S, Jactel H, Jaroszewicz B, Joly FX, Zucker T, Koricheva J, Milligan H, Müller S, Muys B, Nguyen D, Pollastrini M, Raulund-Rasmussen K, Selvi F, Stenlid J, Valladares F, Vesterdal L, Zielinski D, Fischer M (2016) Jack-of-all-trades effects drive biodiversity–ecosystem multifunctionality relationships in European forests. Nat Commun 7(1):11109. https://doi.org/10.1038/ncomms11109

Vidal JP, Martin E, Franchistéguy L, Baillon M, Soubeyroux JM (2010) A 50-year high-resolution atmospheric reanalysis over France with the Safran system. Int J Climatol 30(11):1627–1644. https://doi.org/10.1002/joc.2003

Vilà M, Vayreda J, Comas L, Ibáñez JJ, Mata T, Obón B (2007) Species richness and wood production: a positive association in Mediterranean forests. Ecol Lett 10(3):241–250. https://doi.org/10.1111/j.1461-0248.2007.01016.x

Zeller L, Pretzsch H (2019) Effect of forest structure on stand productivity in Central European forests depends on developmental stage and tree species diversity. For Ecol Manag 434:193–204. https://doi.org/10.1016/j.foreco.2018.12.024

Zhang Y, Chen HYH, Reich PB (2012) Forest productivity increases with evenness, species richness and trait variation: a global meta-analysis. J Ecol 100(3):742–749. https://doi.org/10.1111/j.1365-2745.2011.01944.x

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.