# PHOREAU v1.0: a new process-based model to predict forest functioning, from tree ecophysiology to forest dynamics and biogeography

Tanguy Postic, François de Coligny, Isabelle Chuine, Louis Devresse, Daniel Berveiller, Hervé Cochard, Matthias Cuntz, Nicolas Delpierre, Émilie Joetzjer, Jean-Marc Limousin, et al.

► To cite this version:

Tanguy Postic, François de Coligny, Isabelle Chuine, Louis Devresse, Daniel Berveiller, et al. PHOREAU v1.0: a new process-based model to predict forest functioning, from tree ecophysiology to forest dynamics and biogeography. Geoscientific Model Development, 2025, 18 (20), pp.7603-7679. 10.5194/gmd-18-7603-2025. hal-05121280

HAL Id: hal-05121280

https://hal.inrae.fr/hal-05121280v1

Submitted on 19 Jun 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>

# PHOREAU v1.0: a new process-based model to predict forest functioning, from tree ecophysiology to forest dynamics and biogeography

Tanguy Postic<fnref n="1" />, François de Coligny<fnref n="2" />, Isabelle Chuine<fnref n="1" />, Louis Devresse<fnref n="1" />, Daniel Berveiller<fnref n="3" />, Hervé Cochard<fnref n="4" />, Matthias Cuntz<fnref n="5" />, Nicolas Delpierre<fnref n="3" />, Emilie Joetzjer<fnref n="5" />, Jean-Marc Limousin<fnref n="1" />, Jean-Marc Ourcival<fnref n="1" />, François Pimont<fnref n="6" />, Julien Ruffault<fnref n="6" />, Guillaume Simioni<fnref n="6" />, Nicolas K. Martin-StPaul<fnref n="6" />*, Xavier Morin<fnref n="1" />* (*co-last authors)

Correspondence to: Tanguy Postic (tanguy.postic@cefe.cnrs.fr)

**Abstract.** Climate change impacts forest functioning and dynamics, but large uncertainties remain regarding the interactions between species composition, demographic processes and environmental drivers. While the effects of changing climates on individual plant processes are well studied, few tools dynamically integrate them, which precludes accurate projections and recommendations for long-term sustainable forest management. Forest gap models present a balance between complexity and generality and are widely used in predictive forest ecology, but their lack of explicit representation of some of the processes most sensitive to climate changes, like plant phenology and water use, puts into question the relevance of their predictions. Therefore, integrating trait- and process-based representations of climate-sensitive processes is key to improving predictions of forest dynamics under climate change.

In this study, we describe the PHOREAU model, a new semi-empirical forest dynamic model resulting from the coupling of a gap model (FORCEPS), with two process-based models: a phenology-based species distribution model (PHENOFIT) and a plant hydraulics model (SurEAU), each parametrized for the main European species. The performance of the resulting PHOREAU model was then evaluated over many processes, metrics and time-scales, from the ecophysiology of individuals to the biogeography of species.

PHOREAU reliably predicted fine hydraulic processes at both the forest and stand scale for a variety of species and forest types. This, alongside an improved capacity to predict stand leaf areas from inventories, resulted in better annual growth compared to FORCEPS, and a strong ability to predict potential community compositions.

By integrating recent advancements in plant hydraulic, phenology, and competition for light and water into a dynamic, individual-based framework, the PHOREAU model, developed on the Capsis platform, can be used to

<footnote n="1">¹CEFE, CNRS, Univ. Montpellier, EPHE, IRD, Montpellier, France.</footnote>

<footnote n="2">²AMAP UMR931, Botany and Computational Plant Architecture, INRAE, TA A-51/PS2, Boulevard de la Lironde, 34398 Montpellier Cedex 5, France.</footnote>

<footnote n="3">³Université Paris-Saclay, CNRS, AgroParisTech, Ecologie Société Evolution, 91190, Gif-sur-Yvette, France.</footnote>

<footnote n="4">⁴INRA, UMR 547 PIAF, F-63100 Clermont-Ferrand, France.</footnote>

<footnote n="5">⁵Université de Lorrain, AgroParisTech, INRAE, UMR Silva, 54000 Nancy, France.</footnote>

<footnote n="6">⁶INRAE, URFM, Domaine Saint Paul, INRAE Centre de recherche PACA, 228 route de l'Aérodrome, CS 40509, Domaine Saint-Paul, Site Agroparc, France.</footnote>
<figure><img src="image_4.png" /></figure>

understand complex emergent properties and trade-offs linked to diversity-effects effects under extreme climatic events, with implications for sustainable forest management strategies.

# 1 Introduction

Forests cover approximately 30% of the Earth's land surface, hosting the majority of terrestrial biodiversity. They are crucial carbon sinks (Pan et al., 2011), play a vital role in climate regulation (Chapin III et al., 2008), and provide essential ecosystem services to humans (Nadrowski, Wirth and Scherer-Lorenzen, 2010). However, climate change poses significant risks to forests, including disruptions to forest dynamics (McDowell et al., 2020a), as increasingly extreme environmental conditions have profound effects on forest structure and composition as well as on forest functioning, including massive mortality events (Allen, Macalady, Chenchouni, Bachelet, McDowell, Vennetier, Kitzberger, Rigling, Breshears, E.H. (Ted) Hogg, et al., 2010). Such impacts are assessed through experimental (Gavinet, J. Ourcival and Limousin, 2019; Decarsin et al., 2024) and empirical (McDowell et al., 2020a) studies. Yet, although such approaches are key to understanding and anticipating forests' response to climate change, they cannot cover the entire spectrum of environmental contexts, species compositions, and forest history. By filling those gaps in knowledge, forest models represent key complementary tools to effectively investigate the combined impacts of species composition and climate change on forest dynamics and functioning (Bugmann, 2001; Maréchaux et al., 2021).

Yet the robustness of such models — most often calibrated on historical data — is often questioned when used to make predictions for the uncertain transition period of the coming decades (Parmesan, Morecroft and Trisurat, 2022; Van der Meersch et al., 2025). Focusing on Europe, climate projections generally describe drier conditions, with might lead to a shift from light to soil water as the main limiting resource over which individual trees compete (McDowell et al., 2020a). In this context, the accuracy of forest projections might depend in large part on whether models are able to account for causal relationships between water stress and stand composition (Brodribb et al., 2020; McDowell et al., 2022; Van der Meersch et al., 2025). For example, instead of postulating general a priori species complementarity effects in resource use, process-based modelling must strive to capture how individual trees harness and compete for light and water in natura.

Furthermore, depicting and understanding the role of diversity in ecosystem functioning has been a key focus of ecological studies for at least two decades (Kinzig, Pacala and Tilman, 2002; Hooper et al., 2012; van der Plas, 2019). In forest ecosystems, the importance of the role of diversity — both structural and compositional — on productivity and wood biomass has been firmly established by numerous studies over a wide range of conditions and methods (Nadrowski, Wirth and Scherer-Lorenzen, 2010; Morin, 2011; Paquette and Messier, 2011; Liang et al., 2016; Ratcliffe et al., 2017). In addition, there is some evidence that tree diversity could modulate the resistance and recovery of forest productivity under stress or disturbance (Ammer, 2019; Jourdan, Lebourgeois and Morin, 2019; Schnabel et al., 2021; Blondeel et al., 2024), although the level of consensus varies with the type of stress or disturbance considered (Decarsin et al., 2024; Messier et al., 2022). Yet despite these patterns, there remains a scarcity of data regarding the actual differences in functioning of monospecific and mixed forests, and their relative response to changing climate conditions. In fact, while the diversity-productivity relationship is well evidenced — a global meta-analysis has shown mixed-species stands were on average 25% more productive than their respective
<figure><img src="image_5.png" /></figure>

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

species' monocultures (Zhang, Chen and Reich, 2012) —, data regarding the link between species diversity and the ability to withstand extreme climatic events is more scarce and contradictory. Where some studies have linked forest diversity to a lessened sensitivity of tree growth to drought (Lebourgeois et al., 2013; Anderegg et al., 2018; Serrano-León et al., 2024), others have found this relationship to be strongly context-dependent (Grossiord et al., 2014; Forrester et al., 2016; Jactel et al., 2017), and restricted to dry environments. Moreover, with the rapid shift in climatic conditions, it would be a mistake to assume that the same patterns of diversity-productivity and diversity-resilience relationships used to support the stress-gradient hypothesis (Bertness and Callaway, 1994) will apply in the next decades to newly drought-prone sites, where water resource limitation has not had the chance to shape the co-evolution of the local species over the past millennia. In fact, the same structural and specific complementarities that are currently responsible for increasing the productivity of existing mixed temperate forests through a better usage of the light resource could become a source of vulnerability, as competition for water intensifies proportionally to the density and foliage areas of the stands (Jucker et al., 2014; Haberstroh and Werner, 2022; Decarsin et al., 2024; Moreno et al., 2024).

For these reasons, and because experimenting composition effects in mature forests is especially difficult, the evaluation of diversity effects in forest ecosystems has also increasingly relied on forest models, particularly those based on processes (Bohn and Huth, 2017; Maréchaux and Chave, 2017; Jonard et al., 2020; Morin et al., 2021). Indeed, the prospective power of such models make them key tools in testing various hypotheses on the diversity-functioning link (Maréchaux et al., 2021; De Cáceres et al., 2023a), but also in evaluating forest management practices that incorporate species mixing (Jourdan et al., 2021) and more generally in simulating forest-response to the long-term impacts of climate change (Reyer, 2015).

To improve our ability to forecast the impact of climate change on forests and to better test adaptation solutions related to composition and management, we have thus identified two main shortcomings in forest models: the representation of the hydraulic functioning of trees, and of the interspecific interactions — especially competitive ones.

In fact, there is a lack of knowledge regarding the effects of species mixing on forest resistance and resilience to drought, although trait-data describing the hydraulic functioning of tree species has been steadily accumulating in the last years. A great variety of water-stress adaptation and drought response strategies among species have been identified (Martin-StPaul, Delzon and Cochard, 2017; Choat et al., 2018): these include traits linked to the allocation between transpiring and conducting surfaces, stomatal control and conductance (Johnson et al., 2012), water storage, root-to-shoot ratio, specific leaf area, safety margins (Martin-StPaul, Delzon and Cochard, 2017), and rooting depths (del Castillo et al., 2016). These traits and their variability ultimately account for many of the plant-to-plant interactions responsible for water-competition reduction and facilitation (De Cáceres et al., 2021; Moreno-de-Las-Heras et al., 2023; Moreno et al., 2024; Mas et al., 2024). However, understanding their net impact in existing forests is complicated by environmental and structural variability among stands, and more generally by the fact that the most common available indicators — growth and mortality — integrate over time many processes that are difficult to unravel. Therefore, although the dynamic and integrative effect of species-mixing on medium-term drought-resilience most directly concerns forest management strategies elaborated today, it is difficult to
<figure><img src="image_7.png" /></figure>

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

formulate *a priori* recommendations. Decoupling the effects of hydraulic trait diversity from forest structure (foliage area, tree density) involves significant methodological difficulties (Forrester and Pretzsch, 2015), and is further complicated by the feedbacks between traits and stand structure (Guillemot and Martin-StPaul, 2024), as trees have been shown to adapt hydraulic to the forest structure (Limousin et al., 2012, 2022; Martin-StPaul et al., 2013; Moreno et al., 2024).

125

Furthermore, even disregarding species diversity, the relationship between forest structure, density and productivity is itself poorly understood: there is no consensus on the link between tree-size heterogeneity and productivity (Pretzsch and Biber, 2010; Bourdier et al., 2016; Dănescu, Albrecht and Bauhus, 2016), and while stand density has been statically correlated with increased growth (Reineke, 1933; Forrester, 2014), it is the overall dynamic interactions between these factors that must be understood (Morin et al., 2025). The prohibitive cost of testing all the factors affecting forest functioning (species diversity, stand structure and density, response to climate and soil conditions, effect of management...) in experimental or observational studies further justifies the use of forest ecosystem models (Pretzsch, Rötzer and Forrester, 2017), which are able to replicate *in silico* the complex plant-to-plant interactions that regulate competition for above- and belowground resources, evaluate potential facilitation and competition reduction processes, and integrate them over time in stand structure dynamics that account for trade-offs between drought-resistance and productivity.

130

135

140

Recent gap models (Maréchaux and Chave, 2017; Morin et al., 2021) by explicitly modelling crown sizes and species shade tolerances, have focused on capturing the processes through which canopy packing and spatial niche partitioning can emerge. However, space is not only the dimension through which plant species partition resources – time is also an important vector of asymmetry through which different species can coexist in by exploiting different niches (Gotelli and Graves, 1996). Relative shifts of even a few days in leaf phenology – either through earlier budding or later senescence – have been shown to have major impacts on plant growth, by allowing otherwise shaded understory plants to receive full sunlight (Jolly, Nemani and Running, 2004). As warming climate conditions advances the phenology of most species, increasing productivity (Park et al., 2016) at the expense of additional vulnerability to spring frosts (Lopez et al., 2008), accurately integrating phenological responses of individual species is an important next step in improving the ability of gap models to represent competition for light.

145

150

In addition, phenological processes (including seed production, leaf dormancy and resistance to frost) have been shown to be major factors in determining species distribution (Chuine, 2010). Indeed, while many studies highlight the role of species diversity in forest functioning, it is important not to lose sight of the fact that the presence of a species in a given forest is itself the result of a complex historical process conditioned both by site conditions and species coexistence mechanisms. By directly integrating trait-based phenology, gap models can therefore more accurately capture this dynamic by making species diversity an emerging factor of the modelling framework.

155

160

Here we present a new process-based model called PHOREAU which combines the strengths of three previously published process-based models: ForCEEPS, PHENOFIT, and SurEau. The model was developed in the frame of the Capsis modeling plaform (Dufour-Kowalski et al., 2012), a modular software platform designed to simulate
<figure><img src="image_9.png" /></figure>

the growth and management of forest stands. The model extends the scope of classic gap models by including a detailed representation of plant water use and competition for the water resource as well as a detailed representation of plant phenology and its impact on reproduction and frost leaf damage. The PHOREAU model thus presents a coupling between recent advances in the process-based modelling of plant water relations under conditions of extreme drought (Cochard et al., 2021a; Ruffault et al., 2022) with state-of-the art phenology (Chuine and Beaubien, 2001) and light competition (Morin et al., 2021) models, in an individual-based gap-model able to consider most types of forest structures (Morin et al., 2025) and forest management (Jourdan et al., 2021). The validity of this approach is underpinned by its reliance on species-specific hydraulic, allometric and phenological traits, grounded in decades of experimental research (Leinonen, 1996; Kattge et al., 2020; Cochard et al., 2021a).

The PHOREAU model has been designed to shed light on some of the many pending issues regarding the effects of species diversity on forest functioning, such as the impact of extreme droughts (Piedallu et al., 2023) or the role of complementarity in leaf phenology on growth in mixed stands (Morin, 2011). More generally, the model offers the opportunity to tackle issues ranging from the physiology of individuals to the biogeography of species. Therefore, our multi-stage validation protocol, presented here, involves daily hydraulic processes, yearly productivity, pluri-annual mortality, and long-term species composition.

165

170

175

180
<figure><img src="image_10.png" /></figure>

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

## 2 Presentation of the model

The PHOREAU model builds on three process-based models, which have been presented in previous publications. For the sake of clarity, we have chosen to summarize only the main processes of each model, and to focus on the integration methodology and the new processes allowed by the coupling. Refer to Fig. 1 for a schematic representation of the PHOREAU model, and to Fig. 3 for a breakdown of the coupling between the ForCEEPS, PHENOFIT and SurEau models which constitute PHOREAU.

<figure><img src="image_12.png" /><figcaption>**Figure 1** | Schematic representation of the PHOREAU model. The principle of the three main demographic processes (growth, mortality, regeneration) and competition for light are inherited from the ForCEEPS forest dynamics model. Tree hydraulics and competition for water and tree foliar phenology come from the coupling with the SurEau and PHENOFIT models, respectively.</figcaption></figure>

### 2.1 ForCEEPS: a forest community gap-model

#### 2.1.1 Description of the ForCEEPS model

In PHOREAU, forest dynamic processes (growth, mortality and recruitment) are all managed by the ForCEEPS model (Morin et al., 2021). ForCEEPS (Forest Community Ecology and Ecosystem Processes) is a gap model that relies on a few ecological assumptions to simulate the dynamics of tree establishment, growth and mortality in independent small patches of land, that are aggregated to derive properties at the forest scale. While the model is not spatially explicit at the patch level, it is individual-based: two trees of the same species and the same age can have different growth rates under the same climate, depending on the specific patch-level biotic constraints of light-competition. Derived from the FORCLIM model (Bugmann, 1996; Didion et al., 2009) the ForCEEPS model
<figure><img src="image_13.png" /></figure>

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

was developed with the aim of simulating forest dynamics under a wide range of environmental conditions while limiting the need for prior calibration, and was designed to be equally able to simulate planted, managed, or natural forests (Morin et al., 2020, 2025).

Tree growth is computed at a yearly time-step in two phases. First maximum diameter increment is calculated using an empirical equation shown in Eq. 1, as a function of trunk diameter at breast height at the start of the year, and a maximum species growth rate $g_s$. $b_s$ and $c_s$ are species specific allometric parameters (respectively derived from $H_{max}$ and $s$), and $H_{max}$ the maximum height reachable by that species. Height is directly linked to diameter following another species-specific allometric parameter.

$$ \Delta D_{opt} = g_s \times \frac{D \times \left(1 - \frac{H}{H_{max,s}}\right)}{2 \cdot H_{max,s} - b_s \times e^{(c_s \cdot D) \times (c_s \cdot D + 2)}} \qquad \text{Eq. 1} $$

Then, realized growth is determined from optimal growth after reduction by a series of growth-reduction factors (bounded between 0 and 1) following a modified geometric mean, as shown in Eq. 30.

Drought, growing degree days, and soil reduction factors range from 0 to 1 are determined by site soil and climatic conditions, and modulated by species-specific parameters. The other factors represent biotic constraints related to light availability. $GR_{light}$ represents the immediate effect of competition for light, and depends on the cumulated leaf area above or at the same level as the considered tree. $GR_{Crown}$ represents the long-term effects of crown size reduction on the capacity of trees to grow and assimilate carbon. In the ForCEEPS framework, trees crowns were represented as downwards-pointing triangles; the ratio of crown height to tree height $c_s$ is adjusted based on a factor that decreases from a species-specific maximum to minimum value as the tree experiences increasing shade (see Fig. 2).

Similarly, tree establishment is regulated by winter temperature, growing degree days, light availability, and stand browsing intensity (Eq. 32). The number of potential seedlings for a given species depends on maximum site density and its shade tolerance parameter, shade intolerant having a greater regeneration potential (Eq. 31). The survival of each potential seedling is controlled by a stochastic process itself regulated by the reduction factors listed above. If selected, the sapling is initialized with a DBH of 1.27 cm. Tree mortality is the combination of a stochastic background process combining stand density and tree longevity, and a growth-related mortality that represents stress-caused tree death linked to biotic and abiotic constraints.

A full description of the ForCEEPS model developed on the Capsis modeling platform (Dufour-Kowalski et al., 2012) that was used as a base for this study can be found in Morin et al., (2021). In the following section, we present new developments that have been included in the ForCEEPS model, before the coupling with SurEau and PHENOFIT.
<figure><img src="image_15.png" /></figure>

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

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

## 2.1.2 Improvements to the ForCEEPS model

In anticipation of the coupling with SurEau and PHENOFIT, a number of modifications were made to the ForCEEPS model, focusing on microclimate, light-dependent height plasticity, and improvements to the light-competition module. This proved necessary when integrating transpiration-driven water fluxes, as stand leaf area is one the main driver of embolism in the SurEau model (Cochard et al., 2021b), and preliminary results indicated a poor capability of the ForCEEPS to reproduce observed leaf area indices from stand inventory, in both relative and absolute terms. These refinements are summarized below and in Fig. 2, with more in-depth descriptions in supplementary information (appendices A, B, C, D, E and F).

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

**Light-dependent height plasticity:** ForCEEPS infers tree height from trunk diameter using fixed allometric relationships, limiting its ability to capture site effects and competition-driven height-diameter variations. In reality, understory trees allocate more growth to height, while trees in low-density stands prioritize diameter growth (Oliver and Larson, 1996), especially in shade-intolerant species (Delagrange et al., 2004). Recognizing this, we have incorporated dynamic height growth in PHOREAU, by adjusting height increments based on competition-driven parameters and species shade tolerance parameter.

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

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

**Crown-length reversion:** The PHOREAU model improves the representation of crown length dynamics by allowing crown ratio reversion when light availability increases, unlike the ForCEEPS model, which only permitted decline. This adjustment accounts for the impact of tree death or removal on neighboring trees, enabling canopy recovery, with a yearly crown ratio increase capped at 5% of the difference between the previous ratio, and the potential crown ratio based on light conditions.

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

**Species-dependent crown shapes:** The PHOREAU model improves crown-shape representation by allowing for a greater range of crown shapes than the default ForCEEPS inverse-cone, including ellipsoidal and conical shapes. This in turn allows for a better representation of inter-specific competition, with complementarities arising from differences in crown structure.

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

**Density-dependent light availability:** PHOREAU maintains ForCEEPS' balance between predictive power and computational efficiency by simplifying light dispersion calculations, using a vertical stratification approach without explicit tree positioning. However, this method reduces light competition to a single leaf area index (LAI) value, overlooking horizontal canopy structure and gaps that influence tree growth. To address this, PHOREAU integrates a clumping factor (Ω) into its light extinction coefficient, capturing variations in foliage aggregation and improving realism (Nilson, 1971; Black et al., 1991; Bréda, Soudan and Bergonzini). This approach reflects observed trends, such as the inverse relationship between LAI and light extinction (Dufrêne and Bréda, 1995), and aligns with methods used in remote sensing (Demarez et al., 2008; Chen et al., 2012; Zhu et al., 2018).

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

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

**Incorporation of Specific Leaf Area (SLA):** ForCEEPS crown size allometric relationships, originally calibrated for a few temperate European species (Burger, 1951; Bugmann, 1996), led to inaccurate predictions when applied to a broader range of species, particularly Mediterranean and understory trees. PHOREAU addresses this by
<figure><img src="image_25.png" /></figure>

285 recalculating tree foliage area using species-specific leaf area (SLA) values, improving the model’s ability to represent interspecific differences in drought resistance.

*Microclimate derived from stand-structure:* Forest canopies buffer climatic conditions in the understory, resulting in cooler, more stable daytime temperatures and warmer nighttime temperatures compared to the canopy. This microclimate effect is especially pronounced in dense, structurally complex canopies (De Frenne *et al.*, 2021), helping young understory trees resist drought despite shallow root systems (Forrester and Bauhus, 2016). Because the PHOREAU model integrates fine-scale hydraulic and phenological mechanisms within a forest-structure gap model, it is able to capture these effects of microclimate on plant functioning. In particular, we integrate microclimatic temperatures and vapor-pressure deficits derived from macroclimate data using a statistical model based on a slope and equilibrium approach presented in Gril *et al.*, (2023) and Gril, Laslier, *et al.*, (2023), incorporating patch characteristics like leaf area index (LAI), maximum tree height, and vertical complexity index (VCI). Hourly microclimate temperatures are then used to calculate vapor pressure deficits for transpiration computations, as well as degree-day accumulation for tree growth and regeneration.

290

295

300
<figure><img src="image_26.png" /></figure>

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

<figure><img src="image_28.png" /><figcaption>**Figure 2** | Presentation of the modifications in the light-competition module between ForCEEPS (Morin et al., 2021) and PHOREAU, with a description of the main changes</figcaption></figure>
<figure><img src="image_29.png" /></figure>

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

## 2.2 SurEau: a plant hydraulics model

The SurEau model (Cochard et al., 2021; Ruffault et al., 2022) is a model of the SPA family (soil-plant-atmosphere, Mencuccini et al., 2019), dedicated to model plant response during extreme drought, which describes water flows in a soil, plant and atmosphere system. It was developed with the idea (1) simulating the water status of plants throughout a complete drying sequence going beyond stomatal closure, including plant desiccation and hydraulic failure (Choat et al., 2018); and (2) of being able to be initialized from accessible environmental data (climate, description of the structure of the forest stand by inventory or remote sensing) and hydraulic "traits" at fine taxonomic grains (species, provenance, etc.) which are increasingly available in global databases (e.g. Martin-StPaul, Delzon and Cochard, 2017; Guillemot et al., 2022). The SurEau model uses daily climate data as inputs, which are then disaggregated into hourly values; among its outputs are the time to full stomatal closure, and the hourly level of cavitation of each organ. There are two published versions of SurEau and their detailed presentation can be found in Cochard et al., (2021) and Ruffault et al., (2022). These two versions differ in the complexity of the hydraulic architecture of the plant and the numerical scheme used to solve the equations of transport (Ruffault et al., 2022).

We describe below in a synthetic manner the main principles of the model, the equations used for the coupling, and its implementation in Phoreau. For the purpose of the coupling, we have recently implemented a highly modular version of SurEau into the Capsis platform using Java object-oriented programming, which includes the main aspects of both previous versions of SurEau. The specific functioning of each compartment is elegantly implemented using object-oriented principles, allowing for modularity and clarity in the model design.

SurEau includes principles of forest water balance such as transpiration, rainfall interception, soil evaporation, rain infiltration into different soil layers, and water drainage into deep reservoirs. The specificity of SurEau is to explicitly represent water transport within the tree through a system of resistance and capacitance (Fig. 3). This hydraulic architecture makes it possible to calculate the water status (water potential and water content) at different levels of the tree and the soil. The tree's organs (e.g., roots, trunk, branches, leaves) are represented by a water compartment separated into a symplasm and an apoplasm. The symplasm corresponds to the water reservoir made up of living tissues (parenchyma, phloem, etc.); it is elastic and can exchange water with the vascular system under the effect of tissue volume variations.

The soil-plant-atmosphere system is modeled through different compartments ("hydraulic cells"), considered as "computational entities" and implemented as classes in Java, which are interconnected and exchange water fluxes through specific functions which model ecophysiological processes. This Capsis version builds on the implementation of generic computational entities that we call *SPH* (Soil-Plant-Hydraulic) compartments, which can be attributed a specific type (soil, symplasm, apoplasm). Each type is defined by specific functions to compute water potential and water quantities. These *SPH* compartments can be connected together to build a tree of any possible complexity. The fluxes between cells are determined with Fick's law by using the water potential gradients between cells and their hydraulic conductances. The water content of each cell is therefore described as the result of inflows and outflows; and the water potential of each cell is calculated with the appropriate formulation according to the nature of these cells (soil, symplasm, apoplasm). For the soil a water retention curve
<figure><img src="image_31.png" /></figure>

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

is used (Van Genuchten, 1980). For the symplasm, the law of pressure-volume curves (Tyree and Hammel, 1972), which expresses the relationship between water content and water potential, is used to describe loading and unloading dynamics. These laws can be parameterized using abundant pressure-volume curve data (Bartlett, Soffoni and Sack, 2012). The effect of cavitation is to alter the hydraulic conductance of the apoplasm, and can lead to hydraulic failure. However, cavitation also releases apoplastic water into the transpiration stream, which can temporarily attenuate the drop in water potential (i.e., water stress). Both phenomenons are irreversible (but see Sect. 3.4.2). The percentage loss of conductance (PLC) through vessel embolism is calculated using the water potential of the organ's apoplasm ($\psi_{Apo}$) and an empirical sigmoid function described by species-specific inflexion and slope parameters ($P_{50}$, slope$_{cav}$) as shown in Eq. 2:

$$PLC = \frac{100}{1 + (e^{-\frac{slope_{cav}}{25}} \times (\psi_{Apo} - P_{50}))} \qquad Eq. 2$$

PLC is a key indicator of the risk of mortality by hydraulic failure, and has been elected a key variable for the coupling with ForCEEPS (see Sect. 3.4.2).

The main fluxes from the plant to the atmosphere are the stomatal and the cuticular transpirations. Cuticular and stomatal transpirations are computed using gas-phase conductance, and the vapor pressure deficit between the organ and the atmosphere. The leaf stomatal and cuticular conductance are connected in parallel to produce the leaf conductance, itself connected in series to other boundary and crown conductances to produce the overall canopy conductance. Leaf cuticular conductance varies with leaf temperature and its photosynthetic activity. Meanwhile, stomatal conductance is calculated as the product of a maximum stomatal conductance without water stress $g_{stom,clim_max}$ (which ranges between species specific parameters $g_{stom_max}$ and $g_{stom_night}$ depending on depends on light, temperature, and CO2 concentration), with a regulation factor $\gamma$ based on plant water status, as shown in Eq. 3.

$$g_{stom} = g_{stom,clim\_max} \times \gamma \qquad Eq. 3$$

In particular $\gamma$ represents the degree of stomatal closure between 0 and 1, computed using leaf symplasm water potential $\psi_{LSym}$ and a sigmoid function described by inflexion and shape parameters $\psi_{gs50}$ and slope$_{gs}$ as shown in Eq. 4 (these parameters are themselves derived from species-specific pressure-volume curve parameters $P_{gs12}$ and $P_{gs88}$ : refer to Ruffault et al., (2022), for more details).

$$\gamma = 1 - \frac{1}{1 + e^{-\frac{slope_{gs}}{25}} \times (\psi_{LSym} - \psi_{gs50})} \qquad Eq. 4$$

Numerical resolution of the plant water balance is based either on the explicit or the faster semi-implicit method presented in Ruffault et al., 2022. This first version of PHOREAU v1.0 uses the same simplified tree hydraulic architecture as in Ruffault et al., (2022) and uses the faster and generic semi-implicit solver. Before performing
<figure><img src="image_33.png" /></figure>

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

the coupling, we verified this new implementation could provide exactly the same results as the previous version
under the same initial conditions.

## 2.3 PHENOFIT: a phenology-based distribution model

The PHENOFIT model (Chuine & Beaubien 2001) is a process-based species distribution model for temperate trees which calculates the probability of presence over several years of a given species for a particular set of environmental conditions. This probability is derived from the estimated fitness of an average adult individual of that species, which is itself the product of the probability to survive until the next reproductive season, and the probability to produce viable seeds by the end of the annual cycle. The model assumes that survival and reproduction depend on the synchronization of tree development to seasonal climatic variations, with the plasticity of key phenological events such as leaf unfolding, flowering, fruit maturation, and leaf senescence. The model uses soil data and daily meteorological data (minimum and maximum temperature, rainfall, relative humidity, global radiation, and wind speed) as inputs. It is composed of several sub-models: phenological models for leafing, flowering, fruiting and leaf senescence (for reviews refer to Chuine and Régnière, 2017, and Chuine et al., 2024); a frost injury model (Leinonen, 1996); a survival model; and a reproductive success model calculated as the proportion of uninjured fruits that reach maturation considering photosynthetic ability and the proportion of leaves not killed by frost (Chuine and Beaubien, 2001). A visual representation of the model can be found in Fig. 3.

In PHENOFIT, both the leafing and the flowering dates ($t_f$) are calculated with a two-phase phenology model. In the first phase of endodormancy (Eq. 5), the bud must be exposed to a certain amount ($C_c$) of chilling units ($R_{c,t}$) from the onset of dormancy ($t_0$) in order to break this endodormancy at date $t_1$. In the second phase of ecodormancy, or quiescence (Eq. 7), the bud cells elongate in response to forcing temperatures. They must accumulate forcing units ($R_{f,t}$) until a threshold value ($F_c$) is reached, that corresponds to the leafing or flowering date. The type of response functions to temperature are identical for leafing and flowering, only the parameters of these functions differ between the two. Calculations are done at daily time-step, using mean daily temperatures ($T_d$) and species-specific parameters ($a, b, c, d, e$) as shown in Eq. 5 and Eq. 7. Leaf senescence dates $t_c$ are calculated following the model of Delpierre et al., (2009).

Flowering and leafing dates are then used, alongside the daily minimum temperature ($T_m$) between bud onset and leaf senescence or fruit maturation, to determine proportions of leaves and flower-fruits ($I_t, I_f$) uninjured by frost. The probability that fruits reach maturation ($I_r$) is calculated on the basis of the proportion of uninjured leaves which produce the assimilates accumulated in the fruits, the date of flowering from which thermal energy can begin to be accumulated, and a species-specific parameter $E_c$ representing the average amount of energy needed to reach maturation (Eq. 11). Finally, a yearly probability of producing viable seeds, or reproductive success ($R$), is calculated as the product of the probability that fruits will ripen and the proportion of uninjured fruits reaching maturation, as shown in Eq. 12.
<figure><img src="image_35.png" /></figure>

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

$$C_c = \sum_{t_0}^{t_1} R_{c,t} \qquad \text{Eq. 5}$$

425

$$R_{c,t} = \frac{1}{1+e^{c(T_t-e)^2+d(T_t-e)}} \qquad \text{Eq. 6}$$

$$F_c = \sum_{t_1}^{t_f} R_{f,t} \qquad \text{Eq. 7}$$

430

$$R_{f,t} = \frac{1}{1+e^{a(T_t-b)}} \qquad \text{Eq. 8}$$

$$I_l = f(t_{leafing}, T_i) \qquad \text{Eq. 9}$$

435

$$I_f = f(t_{flowering}, T_i) \qquad \text{Eq. 10}$$

$$I_r = f(t_{flowering}, I_l, E_c) \qquad \text{Eq. 11}$$

$$R = I_f I_r \qquad \text{Eq. 12}$$

440

For each organ and each species, parameters are inferred statistically using time series of phenological observations from native populations (dates of leaf unfolding, senescence, flowering, and fruit maturation) for different sites and different years, or from experimental results found in the literature (resistance of plant organs to frost).

445

As the model simulates one average individual, it does not take into account demography or biotic interactions with other species. It also does not represent the impacts of plant growth on survival and resource allocation, but takes into account the effect of a reduction of leaf area on survival. While it can (by calibrating parameters from phenological data of different provenances) represent the way phenological plasticity can vary from one site to another due to genetic differentiation and eventual local adaptation, we have chosen here to use only one calibration set per species: in other words, we account for the plastic response of a species to varying climate conditions, but not for the genetic differentiation of this response. As a result, species performance may be underestimated at the limits of its distribution due to non-representative parameter estimates.

450

The version of the model used in the study, as well as each species' response parameters, are distributed on the Capsis platform (Dufour-Kowalski et al., 2012). A description of the model can be found in Chuine and Beaubien, (2001).

455
<figure><img src="image_37.png" /></figure>

## 2.4 PHOREAU: the coupled model

### 2.4.1 Model-coupling framework

At the heart of the PHOREAU model is the integration of the ForCEEPS, SurEau and PHENOFIT models. This integration was made possible by the presence of all three models on the Capsis Java platform (Dufour-Kowalski et al., 2012). Two major considerations guided the coupling of the models: avoiding overlapping processes, and minimizing the increase in computing time that might arise when integrating models operating at different time-scales.

In its simplest state, the connection between the three models can be described as follows. Independent PHENOFIT simulations are first run for each species and climate year, whose outputs (dates of leaf unfolding and senescence, probability of reproduction) are then read and fed into the main PHOREAU simulation.

At the beginning of each PHOREAU simulation year, all the trees currently present in the plot are used to initialize a separate SurEau simulation. This simulation lasts exactly one year, using the same daily climate as the main simulation, albeit with a further hourly disaggregation required by the Sureau numerical scheme. In addition to species hydraulic traits parameters (see Ruffault et al., 2022), morphological (i.e. size dependent) variables (including tree volume computed from height and diameter, as well as leaf area, PLC, and light availability), are retrieved directly from the main ForCEEPS simulation; leafing and senescence dates are obtained from PHENOFIT; and the initial state of the soil is retrieved from its state at the end of the previous SurEau simulation for year n-1.

Once a SurEau simulation has been initialized, it proceeds to run for one year. Throughout the simulation data is collected and sent back to the main ForCEEPS simulation to determine the effects of drought stress on growth, mortality, and defoliation, as detailed in the following sections.

However, the sub-hourly time-scales of the SurEau processes, which represent a roughly tenfold increase in computation time, warranted the implementation of two major optional simplifications to this framework. They are summarized below, with more in-depth descriptions in supplementary information (Appendices G and H).

*Treewise aggregation for SurEau module.* SurEau simulation runtimes are primarily influenced by the number of distinct SPH-compartments, and particularly the number of trees. To optimize runtime, PHOREAU reduces the number of trees simulated by SurEau each year, while maintaining the overall stem volumes and foliage areas at the stand, species and cohort level. This is achieved through an aggregation method that groups trees into a predefined number of classes per species (set to 3 in our model evaluation), preserving structural integrity while simplifying competition for water by reducing the number of trees. Trees are distributed into a configurable number of classes based on trunk diameter, separating for example mature and juvenile trees. As trees grow, they may shift between classes, and some classes may remain empty in certain years. Each class is represented by a single aggregate tree, whose characteristics are determined by summing (volume, foliage area, biomass) or averaging (height, root depth, light availability) the corresponding attributes of the individual trees. At the end of
<figure><img src="image_38.png" /></figure>

each year, aggregated class results are distributed among the trees that make them up, informing yearly growth and mortality equations. This method significantly reduces computational complexity, while maintaining key ecological dynamics in SurEau.

<figure><img src="image_39.png" /><figcaption>**Figure 3** | Detailed representation of the processes included in the SurEau, ForCEEPS and PHENOFIT models. Red circles indicate outputs used for the coupling, and red lines their destination in the ForCEEPS simulation. Original figures are taken from Morin et al. (2021), Chuine and Beaubien (2001), and Ruffault et al. (2022), where parameters details can also be found.</figcaption></figure>

### 2.4.2 Drought-stress integration

PHOREAU accounts for drought impacts on tree growth and mortality thanks to the integration of the SurEau plant hydraulics model. Drought-induced mortality can occur either directly — in response to extreme drought through high level of xylem embolism leading to hydraulic failure — or as a long-term consequence of reduced growth related to consecutive low intensity drought and defoliation. As a result, the model effectively represents the interplay between the short term extreme drought effect of hydraulic failure, and the longer term drought effect carbon starvation (McDowell *et al.*, 2008).

*Drought feedback on growth* in PHOREAU is assessed by using the factor of stomatal aperture $\gamma$ computed by SurEau at the tree level. This replaces the ForCEEPS formulation, where a growth reduction factor $GR_{drought}$ was computed by comparing a drought index ($DrI$) based on a simple monthly water budget with an empirical species-specific drought tolerance index (Bugmann and Solomon, 2000). The factor of stomatal aperture $\gamma$ is computed (Eq. 13) from the leaf water potential on the basis of a sigmoid curve described by two species-specific traits ($P_{gs12}$ the water potential causing 12% stomatal closure, and $P_{gs88}$ the water potential causing 88% stomatal closure, Cochard *et al.*, 2021b, Ruffault *et al.*, 2022). Daily stomatal apertures are then integrated annually, over the vegetation period, to compute the $DrI$ (Eq. 14). Refer to Appendix I for more details.
<figure><img src="image_40.png" /></figure>

$$ \gamma = 1 - \left( 1 + e^{\frac{P_{L,sym} - 0.5 \times (P_{gs12} + P_{gs88})}{0.25 \times (P_{gs12} - P_{gs88})}} \right)^{-1} \qquad \text{Eq. 13} $$

530

$$ DrI = 1 - \frac{1}{n} * \sum_{j=1}^{n} (\gamma_j) \qquad \text{Eq. 14} $$

535

n : days in year ; j : day of year

540

Drought feedback on mortality and defoliation. Two additional drought stress mechanisms derived from the level of embolism were implemented in PHOREAU. First, drought-induced defoliation was computed on a daily basis for each tree by using the percentage of the leaf xylem embolism (Cakpo et al., 2024). The defoliation rate was set proportional to the embolism rate, with a minimal threshold set at 10% (Eq. 20 and 21). The resulting defoliation percentage is applied to the maximum leaf area of the tree for the given day (itself the result of the species crown allometry, reduction of crown length due to competition for light, and the phenological stage of the leaf derived from PHENOFIT) to obtain the effective daily leaf areas used throughout the model, from plant water usage to light competition and rain interception (refer to Sect. 3.4.3 for details and equations). Furthermore, an average yearly defoliation percentage is computed for integration in the GRcrown growth-reductor from crown length which represents leaf-loss impact on carbon assimilation (see Eq. 27 to 29). Finally, the longer-term adaptation between water stress and reduced leaf area is partially captured by the fact cavitation is carried over from year to year, with a specific repair mechanism described below. Refer to Appendix J for more details.

545

550

555

Second the rate of embolism (assessed through the percent loss of cavitation, PLC) is used to estimate extreme drought induced mortality. The PLC computed by SurEau is retrieved for each tree at the end of the year. Because no cavitation-repair mechanism is implemented at this intra-yearly timescale, the end-of-year value is also necessarily the maximal reached PLC. Then, the resulting PLC% is converted into a probability of death, which is applied at the end of the year like the other death probabilities in the model (Eq. 33). When the tree aggregation option (see Appendix G) is used, each individual tree of a class receives the drought-induced death probability of its corresponding aggregate tree, and death events are drawn independently among them. The actual conversion of the level of cavitation into a death-probability follows a logistic distribution fitted using data from Hammond et al., 2019. The probability distribution is parametrized using a constant steepness parameter, and a species-specific LD50 parameter which corresponds of a point of no return, the lethal dose of cavitation at which exactly 50% of individuals of the species are expected to die (see Eq. 15). As a first approach this LD50 was fixed parameterized at respectively 50% and 80% for gymnosperm and angiosperm species (Choat et al., 2012); Delzon and Cochard, 2014), reflecting the capacity of the latter species to operate at water potentials below the P50 line. This is a result of differences in strategies between embolism-tolerant and embolism-avoidant species, as gymnosperms tend to operate at wider safety margins with vessels more resistant to embolism (Choat et al., 2012). Finally, an additional threshold parameter was added to avoid random mortality events for low PLC values,

560

565
<figure><img src="image_41.png" /></figure>

considering even well-watered trees show some degree of embolism throughout the year (Cruiziat, Cochard and Amiglio, 2002). Refer to Appendix K for more details.

$$ P_{cavitationMortality} = \begin{cases} (1 + e^{-\lambda*(PLC_{\%}-LD_{50,s})})^{-1} & \text{PLC}_{\%} > PLC_{threshold} \\ 0 & \text{PLC}_{\%} \le PLC_{threshold} \end{cases} \quad \text{Eq. 15} $$

s : species ; PLC%: end-of-year loss of conductance percentage ; LD50,s : species cavitation sensitivity parameter ; λ : steepness parameter (default 0.12) ; PLC_threshold : default 20%

Year-to-year cavitation memory and repair. The impact of cavitation on tree functioning can continue long after the end of the initial drought event, and is one of the main causes for the increased vulnerability to future drought events of previously weakened trees (Anderegg et al., 2013; Feng et al., 2021). On the other hand, internal repair mechanisms linked to plant growth (formation of new vessels) can allow the recovery of initial conductance over time (Brodribb et al., 2010). As such, the recovery from embolism in PHOREAU is driven by basal area growth — or, more precisely, by the relative increase of sapwood area, which contains the living conductive vessels. While all new growth is naturally sapwood, as a tree becomes larger the relative proportion of sapwood to heartwood tends to decreases. It follows that to evaluate the rate of replacement of the conductive vessels, the model must first know the pre-existing area of sapwood. PHOREAU uses the foliage area to determine this quantity, through the application of a species-specific, constant, leaf-to-sapwood ratio, also known as the inverse of the Huber value (Cruiziat, Cochard and Amiglio, 2002). The leaf-to-sapwood ratio is applied to the potential one-sided leaf area of the tree, derived solely from its DBH and allometry parameters, and not its actual leaf area after defoliation through competition, frost or drought. This approach, presented in Eq. 16, assumes the Huber value to be constant: we know that this is in fact an important simplification, and that many species adapt their leaf mass per area to site conditions (Lopez et al., 2008).

$$ PLC_s^{n+1} = Max(0, PLC_s^n - 100 * \frac{\Delta BasalArea_{n+1}}{LAp_{n+1} * LA:SA_s}) \quad \text{Eq. 16} $$

s : species; n : year ; PLC: end-of-year loss of conductance percentage LAp : potential one-sided leaf area ; LA: SA_s: species leaf area to sapwood ratio

### 2.4.3 Leveraging leaf phenology and hydraulics to temporalize competition for light

In ForCEEPS, the way the light availability of each canopy layer is determined by the above total leaf area of the above layers, combined with differentiated shade tolerances between species, allows emergent complementarities in a multi-specific context between shade tolerant and intolerant species, resulting on average in greater total stand leaf area and productivity at the stand level (Morin et al., 2025). But alongside spatial complementarities, there
<figure><img src="image_42.png" /></figure>

exist temporal complementarities in species usage of light related to different leaf phenology (Gotelli and Graves, 1996).

605

The PHOREAU model, by integrating leaf phenology simulated by the PHENOFIT model (see Sect. 3.3), accounts for these temporal effects. In particular, the PHENOFIT model calculates two dates based on temperature and photoperiod conditions: the unfolding date ($t_{f,s,n}$) when 50% of the buds show at least one unfolded leaf (BBCH 15), and the senescence date ($t_{c,s,n}$) when 50% of the leaves have changed color or have fallen (BBCH 95). This gives us the range of days when each tree bears leaves. In practice, the maximum daily foliage area of a given tree ($LA_{s,i}^{n,j}$) is derived from its maximum yearly foliage $LA_{s,i}^n$ (itself the result of species-specific crown allometry and the light availability of the tree, Eq. 17 and 18), by using the dates of leaf unfolding $t_{f,s,n}$ and leaf senescence $t_{c,s,n}$ calculated by PHENOFIT for a given species $s$ for a given year $n$, as described in Eq. 19.

610

Using this information required an in-depth reworking of the light-competition module: instead of calculating each layer's light availability at the yearly time-step, daily light availability is now calculated by summing the crown areas of all leaf-bearing trees in the above layers. The final tree light availability is calculated by summing, over all its layers, for all the days for which it is itself bearing leaves, each daily layer light availability. To correct for the fact that tree growth is dependent on heat as well as sunlight, this sum is weighed using daily mean temperatures. In addition to being temporalized, this formulation integrates all the refinements to canopy representation described in Sect. 3.1.2.

615

620

Furthermore, while ForCEEPS implements a mechanism for competition-driven loss of foliage area, representing the reduction of the crown height of dominated trees as their lower branches die off, it does not incorporate mechanisms of leaf-loss driven by extreme meteorological or hydraulic conditions. Unlike competition-driven branch dieback, leaf-loss caused by extreme weather conditions is not usually accompanied by branch death, does not preferentially target the leaves located in the lower parts of the crown, and can be more quickly reverted with shoot regrowth. These differences justified the implementation in PHOREAU of a new mechanism for transitory leaf-loss, distinct from the reduction of crown size, with no memory from one year to the next. The variables used to drive this leaf-loss are derived from the yearly percentage of frost-damaged leaves ($I_t$) and daily leaf cavitation ($PLC_t$) values calculated respectively in the PHENOFIT and SurEau models (see Sect. 3.2 and 3.3). The PHENOFIT leaf loss index is calculated using the frost injury model of Leinonen (1996), based on the leaf-phenology, temperature and photoperiod conditions. The SurEau drought-induced leaf-loss is presented in Sect. 3.4.4. This new mechanism, shown in Eq. 20 and 21, allows the model to reflect strategies of drought acclimation, where defoliation can help some species tolerate drought events (Bréda et al., 2006; Limousin et al., 2022) at the cost of a lowered growth potential. It is this daily leaf area $LA_{s,i}^{n,j,effective}$ that is *in fine* used in all PHOREAU processes, from transpiration, GDD accumulation for growth, to light-competition. Refer to Eq. 28 and Eq. 20 for the respective formulations of frost-induced leaf loss ($frostComponent_{s,n}$) and drought-induced leaf loss ($droughtComponent_{i,s,n}$).

625

630

635

640
<figure><img src="image_43.png" /></figure>

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

$$LA p_{s,i}^{n} = c_{2}^{s} \times crownsize_{i,n} \times DBH_{i,n}^{A_{2}} \qquad Eq. 17$$

$$ \left\{ \begin{array}{l} c_2^s = 0.35 * SLA_s * 2 \quad (\text{Deciduous}) \\ c_2^s = 0.45 * SLA_s * 2 \quad (\text{Evergreen}) \end{array} \right. $$

$$ crownsize_{i,n} = f(LightAvailability_i) \qquad Eq. 18 $$

$$ LA_{s,i}^{n,j} = \begin{cases} 0 & j \le t_{f,s,n} \\ LA_{s,i}^{n} * \frac{(j-t_{f,s,n})}{UI_s} & t_{f,s,n} < j < (t_{f,s,n} + UI_s) \\ LA_{s,i}^{n} & (t_{f,s,n} + UI_s) \le j \le t_{c,s,n} \\ LA_{s,i}^{n} * \frac{(j-t_{c,s,n}-CI_s)}{CI_s} & t_{c,s,n} < j < (t_{c,s,n} + CI_s) \\ 0 & j \ge (t_{c,s,n} + CI_s) \end{cases} \qquad Eq. 19 $$

$$ droughtComponent_{i,s,n,j} = \begin{cases} 1 - \frac{1}{(t_c - t_f)} \times PLC_{i,s,n,j} & PLC_{\%} > 10\% \\ 1 & PLC_{\%} \le 10\% \end{cases} \qquad Eq. 20 $$

$$ DefoliationPercentage_{i,n,s,j} = frostComponent_{s,n} \times droughtComponent_{i,s,n,j} \qquad Eq. 21 $$

$$ LA p_{s,i}^{n,j,effective} = LA_{s,i}^{n,j} \times \left( 1 - \frac{DefoliationPercentage_{i,n,s,j}}{100} \right) \qquad Eq. 22 $$

s : species; i : tree; n : year; j : day of year; LAp<sub>s,i</sub><sup>n</sup>: maximum tree yearly leaf area; t<sub>f,s,n</sub>: species leaf unfolding date; t<sub>c,s,n</sub>: species leaf senescence date; UI<sub>s</sub>: species leaf unfolding interval; CI<sub>s</sub>: species leaf coloration interval

This simplified formulation has the disadvantage of disregarding intra-specific differences in phenology arising from differences in size or competition-status (Gill, Amthor and Bormann, 1998; Augspurger and Bartlett, 2003; Vitasse, 2013; Gressler et al., 2015). Furthermore, it does not yet take full advantage of the PHOREAU hydraulic submodule to account for the effects of drought on leaf development, either through earlier leaf coloration (Xie et al., 2018) or shifted unfolding (Cleland et al., 2007). Further developments of the PHOREAU model should therefore strive to use information from the light competition and water stress modules to inform the calculation of phenology dates.

The rain interception module. In addition, PHOREAU integrates a rain interception module that reduces incoming rainfall based on daily foliage area, accounting for allometry, competition, frost, phenology, and drought-defoliation effects. Canopy storage volume, derived from daily foliage area, accumulates rainfall and releases
<figure><img src="image_45.png" /></figure>

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

water through evaporation, with throughfall calculated using a simplified Beer-Lambert formula. Refer to Appendix L for more details and model equation.

### 2.4.4 Rooting system representation in PHOREAU

The explicit representation of root and their related processes is crucial for any model aiming to simulate the response of vegetation to climate change (Woodward and Osborne, 2000). Because of this, the framework for representing roots in PHOREAU had to be considerably expanded compared to the parent model where the rooting system was reduced to a simple fine root surface. In particular, we built upon the SurEau-Ecos framework by integrating coarse root depth alongside fine root surface, having the roots of different trees share the same soil to compete for water, and implementing plastic responses of root biomass and root depth to drought stress and aboveground growth.

The modelling of the root compartment in PHOREAU is based on the same major hypothesis as that of the canopy and light competition module: an implicit homogenous horizontal distribution of trees, with an explicit vertical stratification. In the same way the aggregated vertical distribution of foliage area entirely determines the light availability of each tree, competition for soil water between trees in PHOREAU is the result of the vertical distribution of their root systems. The underlying hypothesis is that all trees compete for the same water reserves, provided their roots go deep enough; and the user must take care to select a simulation stand area that verifies this constraint, which will itself vary according to the size and rooting structure of the trees present in the stand.

In PHOREAU the rooting system of a tree is split between fine roots and coarse roots: this distinction is essential as the root types have different functional roles and responses to external factors (Pregitzer, 2002). Schematically, fine roots extend horizontally to absorb water in the available soil, while coarse roots explore deeper layers and make them available to fine root exploration. Because in PHOREAU the soil is segmented in a number of layers, this has been translated in the following way: the fine root area of a tree in a determines the conductance between this tree and a given soil layer, while the rooting depth determines which layers the tree has access to, and how its fine root area is distributed within them.

In practice this means that, for a given set of soil parameters, certain trees will be able to extract water from the full soil profile, while others will be restricted to only a fraction (see Fig. W2, extracted from the PHOREAU evaluation on the ICOS sites). This framework is intended to reflect the crucial role of rooting depth in resilience to drought stress (Canadell et al., 1996), as trees with deeper rooting systems are able to make use of relatively untouched water reserves in deeper soil layers. Furthermore, because this is implemented in a forest dynamics model where many trees share the same soil, PHOREAU will be able to use the differential rooting depths to explore the contrasting intra and inter-specific drought responses observed in nature (Johnson et al., 2018).

Rooting depth is a notoriously difficult trait to measure, and involves costly, time-consuming, usually destructive techniques (Maeght, Rewald and Pierret, 2013). While some rooting depth data is available in the literature (Guerrero-Ramírez et al., 2021), its scarcity makes it difficult to disentangle environmental, allometric, and genetic factors; what is driven by aboveground biomass, from what is driven by water availability and groundwater table
<figure><img src="image_47.png" /></figure>

depth (Fan et al., 2017; Freschet et al., 2021; Li et al., 2022). To circumvent this difficulty in obtaining accurate rooting depth traits, we take advantage of the fact PHOREAU does not explicitly represent the position of a tree in the plot and ignores lateral distribution, by using coarse root biomass — an extensively studied trait — as a proxy for rooting depth, thereby implicitly aggregating the lateral and vertical extension of the root system in an integrative rooting extent variable, which is driven by shoot size and site aridity (Tumber-Dávila et al., 2022).

715

Coarse root biomass and fine root biomass in PHOREAU are calculated independently. Fine root area is derived on a 1:1 basis from leaf area. Meanwhile, coarse root biomass is calculated as the product to above-ground biomass with a root-shoot ratio, this root-shoot itself calculated as ratio of realized tree height to maximum species height, positively modulated by the mean of past drought indices (Morin et al., 2021). This formulation, shown in Eq. 23 to 25, follows the conclusions of Ledo et al., 2018 which identifies size and past droughts as the main factors driving root-shoot. These simple equations allow PHOREAU to capture several well-established characteristics of the evolution of coarse and fine root biomass.

720

725

$$R/S_{realised,s}^n = R/S_{min,s} + \frac{1}{2} * \text{AllometryComponent} + \frac{1}{2} * \text{AdaptationComponent} \quad \text{Eq. 23}$$

$$\text{AllometryComponent} = (R/S_{max,s} - R/S_{min,s}) * \left(1 - \left(\frac{\text{Height}^n}{\text{Hmax}_s}\right)\right) \quad \text{Eq. 24}$$

$$\text{AdaptationComponent} = (R/S_{max,s} - R/S_{min,s}) * \sum_{i=(n-10)}^{n} (\text{DroughtIndex}^i) \quad \text{Eq. 25}$$

730

n : simulation year ; s : species

Similarly to leaf shedding, fine root area tends to decrease in response to past drought events (Hartmann, 2011; Brunner et al., 2015). Meanwhile, total root biomass relative to aboveground biomass (the root-shoot ratio) has repeatedly been shown to be positively correlated to past drought events (Mokany, Raison and Prokushkin, 2006), and tree species adapted to more xeric climates have higher root-shoot ratios and deeper roots than those adapted to wetter conditions. These patterns, captured by PHOREAU (Fig. 4), are in accordance with Optimal Resource Partitioning theory (OPT), which predicts trees should increase their absorptive capacity relative to their transpiring surface under short water supply (Coomes and Grubb, 2000; Hertel et al., 2013).

735

740

Another observation captured by deriving root biomass from relative height in PHOREAU is the negative correlation between root-shoot ratio and above-ground biomass (Mokany, Raison and Prokushkin, 2006; Ledo et al., 2018). Because tree height in PHOREAU tends asymptotically towards the species' maximum height following a parabolic curve, as trees become older they allocate proportionally more growth to their diameter than to their height — and to their roots in the new formulation. Following Konôpka et al., 2010, the maximum root-shoot was set to be greater for angiosperms than coniferous trees, who tend to have shallower roots (Schenk and Jackson, 2002) and less variation between juvenile and adult individuals. Another implication of this formulation is that the proportion of fine roots exponentially decreases with total root biomass (Li et al., 2003). -

745
<figure><img src="image_48.png" /></figure>

750

An emergent property of this framework is that for a given magnitude of water stress, a site which has already suffered past drought events will suffer less mortality and growth loss than a previously wet site, because of the rooting depth adaptation mechanism (Fuchs *et al.*, 2020). This type of plastic adjustment is concurrent with spatial variability in tree dieback related to the level of past drought acclimation (Piedallu *et al.*, 2023). Fig. 4 shows an example of this emergent behavior, by comparing simulations with two different climatic trajectories.

755

This integration of root plasticity, coupled with leaf shedding, is an important first step in the modelling of tree adaptation to drought conditions. However, it by no means provides a complete picture of the various strategies used by trees *in natura*. To refine our approach, the relative importance of past drought conditions relative to that of tree allometry in determining total rooting depth could be determined on a species by species basis, instead of a simple angiosperm/coniferous split. Even then, root plasticity is only one among many plastic responses to drought conditions: regulatory responses have been identified in the ectomycorrhizal network, non-structural carbohydrate concentration, differential gene transcription and pathways, increased suberin and lignin formation in roots, and decreased fine-root turnover rate (Bréda *et al.*, 2006; Brunner *et al.*, 2015).

760
<figure><img src="image_49.png" /></figure>

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

**Figure 4** | Simulated distribution and sum of basal area (m²/ha), LAI and fine root area (m²/ha) at different steps of two 1500 year-long simulation over 50 independent 1000m² mixed *F. sylvatica* and *A. alba* stands, starting from bare ground, using edaphic and climatic conditions from a beech-dominated Vosges stand (RENECOFOR HET 88, 1969-2020). Yearly climates were sampled at random from historical data. For both the control (A) and acclimatization (B) scenarios, incoming rainfall was reduced by 66% from the 1250th year onward. For the acclimatization scenario, an additional 33% reduction was applied between the 1000th and 1250th year.
<figure><img src="image_51.png" /></figure>

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

## 2.4.5 Updated core model equations

Our main concern when coupling the PHENOFIT, SurEau and ForCEEPS models was avoiding that some processes shared by the models be taken into account more than once. For example, we could not directly use the global plant fitness output of PHENOFIT, nor its plant survival output, which integrates drought-effects already represented by the SurEau model. In the end, we used four main yearly PHENOFIT outputs: leaf unfolding and senescence dates ($t_f, t_c$), the percentage of uninjured leaves not damaged by frost ($I_i$) and reproductive success ($R$).

Leaf phenology, i.e. leaf unfolding and leaf senescence dates, were used to control plant fluxes (see Sect. 3.2 and 3.4.3) and the period during which growing degree days (GDD) are accumulated for deciduous species. Evergreen species are assumed to accumulate energy throughout the year. As the ForCEEPS framework worked at a monthly time-step, it was necessary to update the model to calculate GDD using daily temperature data. This introduces both inter-species variability in growth, but also intra-species variability between sites and years. This change impacts both growth (through the temperature growth-reduction factor $GR_{gdd}$) and probability of establishment ($P_{GDD}$). See Eq. 26 for the updated calculation of annual GDD sums, including phenology and microclimate, of a tree of species $s$ and average weighted foliage height $h$, with $T_0$ the base temperature ($T_0 = 5.5^{\circ}C$).

$$GDD_h^s = \sum_{j=t_{f,s}}^{t_{c,s}} \max (T_h^j - T_0) \qquad \text{Eq. 26}$$

Climate-induced leaf loss was integrated into ForCEEPS as a modification to the previously existing crown length reduction factor $GR_{CrownLength}$, which represents the impact of leaf density on growth through carbon assimilation. While in the previous ForCEEPS framework trees could only lose leaves through a lack of light availability (the lightComponent, presented in Morin et al. 2021), PHOREAU also captures drought-induced and frost-damage leaf loss, which are integrated in the updated calculation for $GR_{Crown}$ as shown in Eq. 27 to 29. This is a first approach, following Wang, Zhou and Wang, (2021). We are aware this representation is incomplete, and does not account for leaf regrowth, or differential effects according to tree age and size: the absence of an explicit representation of source and sink compartments, and the lack of tree age data to implement an age-differentiated response to leaf loss, was a limiting factor.

$$GR_{crown} = \text{Min}(\text{lightComponent} \times \text{frostComponent} \times \text{droughtComponent}, 1) \qquad \text{Eq. 27}$$

$$\text{frostComponent} = 1 - \frac{(1-I_i)}{0.5} \qquad \text{Eq. 28}$$

$$\text{droughtComponent} = 1 - \frac{1}{(t_c - t_f)} \sum_{t_f}^{t_c} \begin{cases} PLC_t & PLC_{\%} > 10\% \\ 0 & PLC_{\%} \le 10\% \end{cases} \qquad \text{Eq. 29}$$

$$\Delta D = \Delta D_{\text{opt}} \times \sqrt[3]{GR_{\text{light}} \times GR_{\text{gdd}} \times GR_{\text{drought}} \times GR_{\text{soil}}} \times GR_{\text{crown}} \qquad \text{Eq. 30}$$
<figure><img src="image_53.png" /></figure>

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

Phenofit reproductive success, calculated as the product of the proportion of injured flower-fruits and the proportion of fruits that reach maturity, was used to update the different ForCEEPS regeneration modules. In particular, the yearly number of potential seedlings ($n_{PotentialSeedlings,s}$) for a given species, in addition to its light-tolerance parameter $kLa$, now incorporates its yearly phenology-based reproductive success $R_s$ as shown in Eq. 31. Once the number of potential seedlings for a given species has been determined, the probability of establishment of each individual seedling $P_{est,s}$ is unchanged from the ForCEEPS framework (Eq. 32, with details in Morin et al., 2021), but indirectly integrate the refinements in the calculation of phenology and microclimate (through $P_{GDD}$), light availability at soil level (through $P_{LA}$), and soil water balance (through $P_{Dr}$).

$$n_{PotentialSeedlings,s} = 0.006 \times patchSize_{m^2} \times kLa \times R_s \qquad \text{Eq. 31}$$

$$P_{est,s} = P_{T_w} + P_{GDD} + P_{Dr} + P_{Br} + P_{LA} + c_{est} \qquad \text{Eq. 32}$$

Tree mortality was simply updated to reflect the new cavitation mortality mechanism described in Eq. 15. With $P_0$ and $P_g$ respectively the background and growth-related mortality components described in Morin et al., (2021), the chance that a given tree dies on a given year is such that:

$$P_{mort} = P_{cavitationMortality} + (1 - P_{cavitationMortality}) \times \max(P_0, P_g) \qquad \text{Eq. 33}$$

## 2.5 Model calibration and simulation initialization

*Species parameters.* Species parameters were not tuned on the basis of the evaluation datasets, and, for the majority, correspond to traits determined a priori from the literature and experimental results. A full list of the species parameters used in PHOREAU can be found in Table S13, with accompanying descriptions, examples, and data source references.

*Site parameters.* Site climatic and edaphic conditions were constructed using a mix of on-site measurements, and publicly available European datasets (see Sect. 3.1).

*Crown-length Bootstrapping.* To avoid initial oscillations in stand leaf area resulting from year-wise adjustments of tree crown length based on above leaf area, an algorithm, presented in Appendix M, was developed to initialize all tree crown lengths at equilibrium values at the beginning of the simulation.
<figure><img src="image_55.png" /><figcaption>835 3 Model Evaluation</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>**Figure 5** | Proposed framework for PHOREAU validation. In red the evaluation dataset (described in Sect. 4.1), in green the evaluated model outputs.</figcaption></figure>

840

The key novelty of the PHOREAU model is that it is designed to predict a wide range of forest characteristics and ecosystem functioning features, occurring at various scales. Therefore, we evaluated the model across a broad spectrum of outputs, ranging from daily plant physiological measurements to long-term species composition predictions. This comprehensive approach allowed us to avoid one of the common pitfalls of gap-model, which are often validated on a single integrative metric — such as predicted total stand basal area — which limits the robustness of the predictions under future conditions. By directly assessing the model’s ability to reproduce intermediary variables, such as leaf area indices or soil water fluxes, we could control for common biases that may arise from errors offsetting each other under current conditions, which may not hold true when projecting into future climatic scenarios.

845

Depending on the targeted variable (and especially the available data to characterize it), the model evaluation was conducted on certain sites in France, or on many sites over Europe. Because PHOREAU is intended to be continuously improved and refined over time, the validation protocol and all associated data — summarized in Fig. 5 — will serve as a baseline to evaluate any future modifications to the model.

855
<figure><img src="image_57.png" /></figure>

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

860

## 3.1 Data sources

### 3.1.1 ICOS sites

865

We used data from the Integrated Carbon Observation System (ICOS) for our most in-depth validation protocol that includes hydrological, growth, and mortality components. In particular, we selected four forested sites from the terrestrial ICOS Ecosystem network: Puéchabon, Font Blanche, Hesse, and Barbeau. Together these sites represent a diversity of the climatic, edaphic, and biotic conditions that can be found in France (Fig. 7). Refer to Appendix N for general details on the ICOS network.

870

A preliminary task was building an exhaustive database of all relevant input and output variables over the selected sites. This was made possible by the collaboration of each of the site PIs, especially for non-flux data that was not always readily available on the ETC database (Reichstein et al., 2005; Papale et al., 2006). Table 1 provides a summary of the ICOS data sources used in the model evaluation, as well as some of the main site characteristics, while a more in-depth description of each site can be found in supporting information (Appendices O, P, Q and R). Eventual gaps in data were corrected by selecting, for each of our four sites, the simulation period where the most harmonized data was available. Fig. 6 shows a simulated representation of the initial state of each inventory, highlighting the structural diversity across sites, and Fig. W1 a vertical representation of leaf area distribution.

875

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

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

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

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

<figure><img src="image_63.png" /><figcaption>Figure 6 | 3D vizualisation of ICOS stands used for in-depth validation.</figcaption></figure>

Visualisation generated by the PHOREAU model, on the basis of initial inventories.

880
<figure><img src="image_64.png" /></figure>

<table><caption>**Table 1** | Selected stand characteristics for the four ICOS sites used in the in-depth PHOREAU validation, with associated data sources.</caption><thead><tr><th></th><th>Barbeau</th><th>Font-Blanche</th><th>Hesse</th><th>Puéchabon</th></tr></thead><tbody><tr><td>Location</td><td>48°28′N, 2°46′E<sup>1</sup></td><td>43°44'29"N, 3°35'45"E<sup>1</sup></td><td>48°40'30"N, 7°3'59"E<sup>1</sup></td><td>43°14'27"N, 5°40'45"E<sup>1</sup></td></tr><tr><td>Altitude</td><td>100 m above sea level<sup>1</sup></td><td>425 m above sea level<sup>1</sup></td><td>300 m above sea level<sup>1</sup></td><td>270 m above sea level<sup>1</sup></td></tr><tr><td>Simulation Period</td><td>2006-2021</td><td>2007-2020</td><td>1999-2010</td><td>2003-2020</td></tr><tr><td>Simulation Patch Area</td><td>9 x 1000 m<sup>2</sup></td><td>24 x 267 m<sup>2</sup></td><td>4 x 300 m<sup>2</sup></td><td>3x 100 m<sup>2</sup> (MIND control plots)<sup>4</sup></td></tr><tr><td>Stand Inventory</td><td>Basal area aggregated by size and species<sup>1</sup></td><td>Individual DBH measurements<sup>1</sup></td><td>Individual DBH measurements<sup>1</sup></td><td>Individual DBH measurements<sup>1</sup></td></tr><tr><td>Mean annual temperature</td><td>11.2°C<sup>1,3</sup></td><td>14.8°C<sup>1</sup></td><td>10.1°C<sup>1</sup></td><td>13.6°C<sup>1</sup></td></tr><tr><td>Mean annual precipitation</td><td>677 mm<sup>1,3</sup></td><td>703 mm<sup>1</sup></td><td>948 mm<sup>1</sup></td><td>987 mm<sup>1</sup></td></tr><tr><td>Soil Description</td><td>Endostagnic lovisol over calcareous bedrock<sup>4</sup></td><td>Silty clay loam<br>50%-90% rock fraction<br>Limestone bedrock<sup>2</sup></td><td>Luvic cambisol with local stagnic tendencies<br>Deep loam clay layer<sup>1,2</sup></td><td>Silty clay loam<br>75%-90% rock fraction<br>Limestone bedrock<sup>3</sup></td></tr><tr><td>Available Soil Water Quantity (over 5m profile)</td><td>405.3 mm<sup>3</sup> (extrapolated)</td><td>178.4 mm<sup>1</sup></td><td>447.9 mm<sup>2,3,4</sup> (extrapolated)</td><td>130 mm<sup>1</sup></td></tr><tr><td>Dominant tree species</td><td>Sessile Oak (Quercus petraea)<br>European hornbeam (Carpinus Betulus)<sup>1</sup></td><td>Aleppo pine (Pinus halepensis Mill.)<br>Holm oak (Quercus ilex L.)<sup>1</sup></td><td>European beech (Fagus Sylvatica L.)<br>European hornbeam (Carpinus Betulus)<br>Silver birch (Betula Pendula)<sup>1</sup></td><td>Holm oak (Quercus ilex L.)<sup>1</sup></td></tr><tr><td>Initial Basal Area</td><td>25.4 m<sup>2</sup> / ha<sup>3</sup></td><td>19.6 m<sup>2</sup> / ha<sup>3</sup></td><td>19.4 m<sup>2</sup> / ha<sup>3</sup></td><td>30 m<sup>2</sup> / ha<sup>3</sup></td></tr><tr><td>Dominant Tree Height</td><td>25 m<sup>1</sup></td><td>Pine: 13.5 m<sup>1</sup><br>Holm Oak: 5.5 m<sup>1</sup></td><td>18.3 m<sup>1</sup></td><td>5.5 m<sup>1</sup></td></tr><tr><td>Initial Stem Density</td><td>212 / ha<sup>2</sup></td><td>1008 / ha<sup>3</sup></td><td>3297 / ha<sup>1</sup></td><td>4900 / ha<sup>3</sup></td></tr><tr><td>Stand thinnings</td><td>2011: 15% of basal area<sup>1</sup></td><td>No</td><td>2005: 25% of basal area<br>2010: 15% of basal area<sup>1</sup></td><td>No</td></tr><tr><td>Leafarea index (LAI)</td><td>3.5 - 6.4<sup>1,2</sup></td><td>2.9<sup>2</sup></td><td>4.6 - 7.6<sup>1</sup></td><td>2.2<sup>2</sup></td></tr><tr><td>Flux data</td><td>Provided by Site PI</td><td>Provided by Site PI</td><td>Provided by Site PI</td><td>Provided by Site PI</td></tr><tr><td>Tree water potentials</td><td>Provided by Site PI</td><td>Provided by Site PI</td><td>Betsch et al., (2011)<br>Peiffer et al., (2014)</td><td>Provided by Site PI</td></tr><tr><td>References</td><td><sup>1</sup>: Delpierre et al., 2016<br><sup>2</sup>: Briere et al., 2021<br><sup>3</sup>: Maysonnave et al., 2022<br><sup>4</sup>: Site PI</td><td><sup>1</sup>: Monero et al., 2021<br><sup>2</sup>: Simioni, Marie and Huc, 2016<br><sup>3</sup>: Site PI</td><td><sup>1</sup>: Granier et al., 2008<br><sup>2</sup>: Dufrene et al., 2005<br><sup>3</sup>: Granier et al., 2000b<br><sup>4</sup>: Tóth et al., 2017<br><sup>5</sup>: Site PI</td><td><sup>1</sup>: Limousin et al., 2012<br><sup>2</sup>: Limousin et al., 2022<br><sup>3</sup>: Rambal et al., 2014<br><sup>4</sup>: Gavinet, Ourcival and Limousin, 2019<br><sup>5</sup>: Site PI</td></tr></tbody></table>

885
890

### 3.1.2 ICP II sites

To evaluate our model's predictions of tree and stand productivity, potential natural vegetations, and observed foliage areas, we used 250 plots spread across Europe, from 37.03° N to 69.58° N, and 8.17° W to 30.71° E, covering most of the major European species (Fig. 7. They cover a large range of environmental conditions, with mean annual temperatures (MAT) ranging from -1.62 to 17.6 °C, mean annual precipitation sum (MAP) ranging between 405 and 2707 mm, growing degree days (GDD) ranging from 475 to 4287 °C, and available water quantities ranging from 30 to 671 mm over the soil profile. Refer to Fig. 17 for the distribution of site abiotic conditions, and Table S2 for a detailed site by site breakdown.

895

900

The RENECOFOR network. Following the framework of the ForCEEPS validation (Morin et al., 2021), the RENECOFOR permanent forest plot network was used as the primary validation dataset (Ulrich, 1997). RENECOFOR makes up the French portion of the European ICP II network. Comprised of 102 plots (ca. 0.5 ha) in even-aged managed forests, each composed mostly of a single dominant species, they cover most of the main tree species and environmental conditions in France — with the notable exception of Mediterranean conditions. From the year 2000 onwards, the plots were exhaustively inventoried every five years, as well as before and after

905
<figure><img src="image_65.png" /></figure>

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

every eventual thinning. After the removal of the plots that had suffered the strongest perturbations — and in particular the 1999 windstorm — 97 plots remained. With these, we constructed 192 testing datasets, by grouping for each plot between 2000 and 2021 every pair of inventories that were separated by a period of at least four years within which no disturbance was recorded. The mean initial basal area of the plots was 28.3 m²/ha, while the time-interval between inventories ranged from 4 to 15 years, averaging at 7.1 year. As a rule, we avoided longer time-lapses, which would have mechanically improved simulation results, while giving less information on true model performance.

910

915

*The ICP II network.* In addition, we also used 148 plots from the International Co-operative Program on Assessment and Monitoring of Air Pollution Effects on Forests (ICP Forests), which comprises a network of intensively monitored forest sites (level II plots) distributed across Europe (de Vries *et al.*, 2003; Schwärzel et al., 2022). These plots, located in various European countries, allowed the testing of the model over a wider range of abiotic and biotic conditions. This extension of the validation protocol was facilitated by the fact the RENECOFOR network is the French declination of the European-level ICP II program, with comparable protocols and measurements. Unlike for RENECOFOR, each plot corresponds to exactly one simulation dataset, with no repeat inventories separated by intervals of years. The mean initial basal area of the plots was 28.1 m²/ha, while the time-interval between inventories ranged from 2 to 10 years, averaging at 4.6 years (refer to Table S2 for details on each individual simulation dataset).

920

925

<figure><img src="image_67.png" /><figcaption>**Figure 7** | Spatial distribution of sites used for PHOREAU validation. Sites are color-coded based on the dominant species identified in the inventory (see legend in top-left). Red-bordered diamonds represent the four ICOS site (Puéchabon, Font-Blanche, Barbeau, and Hesse) selected for in-depth hydraulic validation.</figcaption></figure>
<figure><img src="image_68.png" /></figure>

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

### 3.1.3 Climate and soil data

PHOREAU requires detailed daily climatic inputs, as well as comprehensive soil moisture retention measurements (see Table 1). To evaluate PHOREAU we used the ERA-5 Land dataset, a climate reanalysis providing various fields over the world at ~9km resolution (Muñoz-Sabater et al., 2021). The hourly data was aggregated to produce daily time-series from 1969 to 2021 over Europe for our study. Potential evapotranspirations were then calculated at the same resolution using the Penman-Monteith equation (Monteith, 1965).

PHOREAU requires, for each layer of soil (in this study 30 layers, up to a total depth of 5m, see Sect. 3.4.4), the fraction of coarse elements, as well as the parameters of the Van Genuchten water retention curve which describes the soil texture (Van Genuchten, 1980). These parameters were obtained for several depths from the European Soil Hydraulic Database (ESDAC) (Tóth et al., 2017), and interpolated over the height of the soil profile.

The resulting ESDAC soil and ERA-5-Land climate parameter files were used as a baseline for our European validation, and were directly used for the ICP II plots, for which no other climatic or soil data was available. When possible, we completed this continental-scale data with higher-resolution measurements. Field measurements were available for all four ICOS sites, as well as for the RENECOFOR plots for which we used a combination of soil measurements and the SILVAE climate time-series to refine our validation. The mean-correction method used to integrate daily ERA-5 and monthly SILVAE climate time-series is presented in Appendix T. The workflow for climate reconstruction is summarized in Fig. 8.

On-site climate measurements were available for 26 of the 102 RENECOFOR sites (see Table S4 for the list of sites). For some of the sites the measurement periods only partially matched the simulation periods, while for others they were continuous from 2000 to 2021. These datasets, although not directly used in our evaluation protocol (so as not to bias our results for certain sites and species) were instead used to validate our climate reconstruction: first through direct comparisons of climate variable means and variances, and then by comparing the outputs of the ForCEEPS simulations carried-out with on-site vs. reconstructed climatic data (refer to Table S13).

Local measurements of SWHC were available up to a depth of 1 meter for all RENECOFOR plots (Brethes and Frankreich, 1997). Additional measurements were available up to 2 meters for more than half of the plots (Brethes and Frankreich, 1997; Lebourgeois, 2006; Guillemot, unpublished data), which were used to refine validation soil parameters.
<figure><img src="image_70.png" /></figure>

<figure><img src="image_71.png" /><figcaption>**Figure 8** | Summary of the workflow used for constructing PHOREAU evaluation inventories and climate datasets.</figcaption></figure>

## 3.2 Evaluation Protocol

### 3.2.1 Evaluation against intra-annual stand fluxes and tree hydraulics

For each of the four ICOS sites, model predictions were compared to observations at two distinct levels. First for stand-level structure, focusing on the annual trends of leaf area, basal area, and tree mortality, for which statistical metrics were not applied, but predictions instead served as a baseline to identify discrepancies between observed and predicted fluxes (but refer to Sect. 4.5 and 4.6 for direct evaluations on stand productivity and leaf area).

Second, for stand fluxes and tree functional dynamics, measured at the daily level. The performance of the PHOREAU model in reproducing the hydraulic functioning of forest stands was assessed for the following variables (from the most aggregative to the most specific): stand real evapotranspiration (ETR); evolution of soil water content (SWC); tree transpiration derived from sapflow; and stem water potential. Model performance was assessed using the Pearson correlation coefficient (r), the root mean square error (RMSE) and the mean deviation (MD) between observations and model predictions.

### 3.2.2 Evaluation against leaf area index

The evaluation of PHOREAU's ability to predict leaf area indices from inventories was realized on two different levels: first, by comparing model results to those obtained from satellite data for 340 sites spread over Europe featuring a large range of tree species; second, by comparing model results to LAI observations inferred from litter retrieval experiments for a few dozen sites in France.
<figure><img src="image_72.png" /></figure>

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

The novelty of this kind of validation, as well as its importance when considering the fact PHOREAU predicts plants water use without any *a priori* fixing of foliage area (unlike most other tree hydraulics models), are presented and discussed in further detail in Appendix U.

The LAI satellite data used was retrieved from the Copernicus Global Land Service time series derived from daily PROBA-V satellite observations between 1999 and 2020 — first at a 1km resolution, then at 300m from 2014 (Fuster et al., 2020). For all RENECOFOR and ICP II sites and dates used for productivity validation (see Table S2) we compared LAI values predicted from the inventories at the start of the simulation, to those observed by PROBA-V and averaged over the summer months of the given year (but note these values are themselves uncertain (Fang et al., 2019) and likely underestimated for the denser sites).

LAI evaluation on litter data was restricted to those RENECOFOR sites where such data was available — mostly beech and oak sites, excluding coniferous-dominated stands not suited to litter retrieval (Ulrich, 1997).

### 3.2.3 Evaluation against productivity

For each of the 340 selected RENECOFOR and ICP II simulation plots, five patches of 1000 m² were initialized using the inventory of the first inventory campaign (see Table S2). For each patch, trees were sampled at random within the first inventory, until the basal area per hectare of the simulated patch matched that of the original inventory. Sampling was done without repetition within each patch, but with repetition among patches. Trees that were absent from the second inventory or found dead were kept in the sampling in order to match simulated plots to real inventories, but were removed after for growth comparison. As the time step for validation was deliberately kept short, model mortality — either due to stress, age or density — were deactivated for this productivity validation protocol, so as to have for each sampled tree the observed and simulated final diameter. To benchmark model performance, PHOREAU simulation results were compared against ForCEEPS predictions.

For tree species currently not parametrized for ForCEEPS (see Table S13 for a list of the 35 parametrized species), such as *Pyrus communis* or *Ilex aquifolium*, we used one of the generic sets of parameters. In addition to mortality, seedling regeneration was also deactivated in the model, due to the short time scales considered. The *crown A1* ratio between tree height and foliage height was initially set at the species maximum value, and initialized with the canopy bootstrap algorithm (see Fig. M1).

Simulations were run for each site over the time periods indicated in Table S2, repeated five times for each of the five sampled patches. We compared simulated and observed basal area growth at both the tree scale and the stand scale, using predicted and observed basal area increments (BAI) normalized to mean annual values. While comparing actual, instead of averaged, annual increments would have constituted a stronger test, this data is not available for size of plots and the range of species and climatic conditions considered here. For stand-level comparisons, results were directly averaged over the five patches. The performance of both the PHOREAU and
<figure><img src="image_74.png" /></figure>

ForCEEPS model were assessed using the Pearson correlation coefficient (r), the root mean square error (RMSE), the average bias (AB), and the average absolute bias (AAB) between observations and model predictions.

### 3.2.4 Evaluation against potential natural vegetation

To evaluate the model’s ability to predict forest composition through long term simulations for a broad range of climatic conditions — thus integrating the effects of all the different processes for mortality, reproduction, phenology, microclimate buffering effect, and competition not directly captured by shorter-term validations protocols —, we compared community compositions simulated by PHOREAU with the predicted potential natural vegetation (PNV) along an environmental gradient. Here, similarly as in Bugmann (1996) and Morin et al. (2021), potential natural vegetation is simply defined as the assumed dominant tree species, assuming no large disturbances, in a space spanned by mean annual precipitations (MAP) and mean annual temperatures (MAT), following Ellenberg (1986), Rameau *et al.* (2008), and San-Miguel-Ayanz *et al.* (2016). For this validation, we used the same 250 sites (RENECOFOR and ICP II) used for the productivity validation, spanning across all the different PNV conditions described in Ellenberg (1986) (Fig. 17).

For each of the 250 sites, we ran 2000-year simulations starting from the bare ground. This simulation length – accounting for seedling establishment, tree growth and mortality – was necessary to ensure the communities were no longer in a transient phase, and had reached the final stage of forest succession with a pseudo-equilibrium composition. The 2000-year climate time series was obtained by randomizing the years for which climatic data was available (1969-2020), which preserved inter-annual variability in climate, but avoided any cyclic trend. For each site we considered 50 independent patches of 1000 m². At the end of the simulation, aggregate species basal areas per hectare were extracted for each simulated site, and compared to assumed PNV dominant species.
<figure><img src="image_75.png" /></figure>

# 4 Results Analysis

## 4.1 Evaluation of water balance and plant hydraulic functioning

The results of the in-depth evaluation of PHOREAU at the four highly instrumented ICOS sites demonstrated a good ability of the model of the model to reproduce observed ecophysiological and dendrological data across a wide range of scales. The model closely followed observed trends in stand basal area (average R² of 0.59, see Table S15), despite the inherent challenge of predicting individual tree mortality (Fig. W3). It accurately captured both the magnitude and variability of dieback across sites, in terms of both tree density (Fig. 9) and basal area loss (Fig. W5), with a marked increase in the rate of basal area loss in the latter years of each simulation; however, the model slightly overestimated mortality numbers on average and particularly at Hesse (+ 56%, see Table S17), as well as the share of large tree death relative to medium trees and saplings. Predicted foliage area results aligned well with observations in the two open evergreen sites with low mean deviations (*Puéchabon*: 0.19; *Font-Blanche*: 0.35, Fig. W4, Table S16). PHOREAU captured the quick regrowth in foliage area observed at Hesse after the 2005 cut (Granier *et al.*, 2008); however, when comparing absolute values, PHOREAU noticeably underestimated foliage area in the two denser deciduous forests, consistent with prior validation results on leaf area (see Sect. 5.3). Despite these biases, the overall alignment between predicted and observed forest dynamics provides a solid foundation for comparing stand functioning and tree physiological responses at fine temporal resolutions.

The PHOREAU model predicted daily evapotranspiration (ETR) across three of the four ICOS sites, with relatively low mean deviations (*Puéchabon*: 0.03; *Barbeau*: -0.24; *Hesse*: 0.8) and good Pearson correlations (*Puéchabon*: 0.64; *Barbeau*: 0.79; *Hesse*: 0.62) between observed and predicted values (Fig. 10 and V6). At *Font-Blanche*, correlation was moderate ($r$ = 0.48, $p$ < 0.001), as the model underestimated summer ETR while overestimating winter and autumn ETR. This discrepancy, particularly the underestimation of *Q. ilex* transpiration (Fig. W8), may stem from biases in the model's repartition of leaf area between *Q. ilex* and *P. halepensis* and a dampened response of *P. halepensis* stem water potential to summer drought (Fig. 12). Over time, across all sites, the differences between predicted and observed monthly cumulative ETR became more pronounced, reflecting a drift between predicted and observed forest structure. The model also underestimated ETR during the leafless winter months at *Barbeau* and *Hesse*, which could result from the exclusion of understory shrubs from the simulations.

PHOREAU consistently demonstrated good performance in predicting the daily evolution of soil water content (SWC), with low mean deviations (*Puéchabon*: 15.4; *Font Blanche*: 1.03; *Barbeau*: -47; *Hesse*: -31.4; Table S12, Fig. W7) and high Pearson correlations (*Puéchabon*: 0.8; *Font Blanche*: 0.86; *Barbeau*: 0.92; *Hesse*: 0.78) between observed and predicted values. The model generally captured the seasonal refilling of soil water reserves well (Fig. 11). However, at *Hesse*, predicted SWC noticeably lagged behind observations: this is consistent with the model's overestimation of *F. sylvatica* water stress during the 2003 drought (Bréda *et al.*, 2006), and the overestimation of mortality and post-2003 stand ETR (Fig. 11). The possible existence of a temporary aquifer present at the site that was not represented in the model may likely contribute to these discrepancies (Joetzjer & Cuntz, pers. comm.).
<figure><img src="image_76.png" /></figure>

The model also accurately captured the variability in measured leaf water potential across species, seasons, and times of day (Fig. W8). It achieved strong correlations between observed and predicted values for both daily minimum stem potential ($r = 0.71, p < 0.001, n = 208$; Table S10) and predawn stem potential ($r = 0.79, p < 0.001, n = 303$; Table S9), with fair levels of prediction accuracy (RMSE = 0.92 and 0.89, respectively). Despite these strong correlations, the model tended to attenuate the range of observed potentials, underestimating predawn potentials (MD = -0.5) while simultaneously overestimating minimum potentials (MD = 0.53). This bias was particularly noticeable in the predawn potentials of *F. sylvatica* (MD = -1.5), likely in link with the lag in the soil water refilling, though the overall strong correlation ($r = 0.99$; Table S9) highlight the model's ability to reproduce relative trends in tree stress.

1090

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

1095

<caption>**Figure 9** | Predicted versus observed annual tree mortality. For each simulation site, the bars depict the total annual number of dead trees, irrespective of cause, broken down by species and size class contributions (refer to Annex X for details). Observed values are derived from stand inventories, while predicted values are generated by the PHOREAU model. Also shown are the annual mortality rates, calculated relative to the initial number of trees for two distinct time periods in each simulation, along with the total number *N* of dead trees by hectare. Transparent bars indicate years with thinnings (see Table S17 for details), which are excluded from the mortality statistics.</caption>

1100

<figure><img src="image_78.png" /></figure>
<figure><img src="image_79.png" /></figure>

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

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

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

<figure><img src="image_83.png" /><figcaption>**Figure 10** | Predicted versus observed evolution of monthly real Evapotranspiration (ETR). For each simulation site, the bars depict the monthly ETR (mm) predictions generated by the PHOREAU model, broken down by source of flux. Soil and intercepted water evaporation respectively originate from the first layer of soil and the water stored on the surface of leaves, while the two other sources are transpiration from different compartments of the PHOREAU tree (refer to Table S11 for details). The black points indicate the observed monthly actual evapotranspiration (with interpolated lines) representing the total water vapor released from the soil and vegetation into the atmosphere, aggregated from hourly or sub-hourly measurements obtained from each site's flux tower. For the Hesse site, observed ETR has been upscaled from measured sap flux densities.</figcaption></figure>

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

1105

1110
<figure><img src="image_85.png" /></figure>

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

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

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

<figure><img src="image_89.png" /><figcaption>**Figure 11 | Predicted versus observed evolution of soil water content (SWC).** For each simulation site, the black points indicate the observed daily actual SWC, with interpolated lines. The stacked bars depict the daily SWC (mm) predictions generated by the PHOREAU model, with individual contributions of each soil layer stacked and color-coded by soil layer (see Fig. W2 for layer details, and Table S12 for statistics). The predictions are confined the maximum measured depth for each site, as indicated in the upper right corner of the figure. For *Barbeau* and *Font Blanche*, observed SWC were directly obtained from site PIs; for *Puéchabon* and *Hesse*, they were interpolated from soil relative humidity (RH%) measured at different depths, using the same rock fractions as used in the simulation.</figcaption></figure>

1115

1120

1125

1130

1135

1140
<figure><img src="image_90.png" /></figure>

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

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

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

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

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

**Figure 12** | Evolution of predicted versus observed stem water potentials. For the dominant species of the four ICOS simulations, the blue line depicts the daily evolution of the stem water potentials (mPa) generated by the PHOREAU model and averaged over the aggregate trees of the species (refer to Appendix G for details on the aggregation method). The red points represent the observed water potentials, limited to the years for which observational data is available (data sources are detailed in Table 1, and associated statistics in Table S10). For Puéchabon, Font Blanche and Barbeau sites, the minimum daily observed and predicted water potentials are shown. For Hesse, where only predawn observations are available, the maximum predicted water potential is used instead.
<figure><img src="image_96.png" /></figure>

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

<figure><img src="image_98.png" /><figcaption>Figure 13 | Aggregated predicted versus observed daily stem water potential.</figcaption></figure>

All available stem water potentials (mPa) observations are plotted against the PHOREAU predictions for the corresponding day and species. For each species, the full colored lines are the regression lines of the linear model of the relationship between observed and predicted minimum water potential, with confidence interval represented with the grey dashed lines. The dashed red line is the 1:1 line. (a) Comparison with minimum water potentials. (b) Comparison with predawn water potentials.

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

1145

## 4.2 Evaluation against leaf area index

Beyond local litter-based measurements, PHOREAU also demonstrated a reasonable capacity to estimate stand leaf area index (LAI) from observed data across many species and site conditions throughout Europe. When compared to PROBA-V satellite data (Fig. 14), the model yielded a good correlation between observed and predicted LAI values ($r = 0.55, p < 0.001, n = 340$; Table S6), with acceptable prediction accuracy (RMSE = 1.41, AB = 0.08). Although no significant systematic bias was detected, the model tended to dampen the observed variability in LAI, slightly underestimating LAI in denser forest canopies while overestimating it in more open plots.

1150

1155

A species-specific analysis revealed notable biases for certain species. The model consistently overestimated the LAI of dense coniferous plantations, particularly for species such as *P. abies*, *A. alba*, and *P. menziesii*. Conversely, it significantly underestimated LAI for low basal area inventories dominated by *P. pinaster*, *P*.
<figure><img src="image_100.png" /></figure>

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

sylvestris, which could partially result from discrepancies between inventories observed and simulated before and after thinnings. Overall, while PHOREAU presents a notable improvement in capturing inter-species LAI variability compared to the ForCEEPS model (RMSE = 3.42, AB = 0.49; Table S3), it proved less effective in predicting small variations in LAI among structurally similar plots dominated by the same species. When comparing predicted LAI to those inferred from litter collections (for a smaller subset of oak and beech-dominated sites where such data was available) the model did not exhibit any significant bias (RMSE = 0.65, AB = -0.03, n = 40; Table S7, Fig. W12), but showed only middling predictive power ($r$ = 0.3, $p$ = 0.047; Table S7). While this evaluation is necessarily hampered by the fact the observed PROBA-V LAI are themselves reconstructed from reflectance values collected at a 300m² scale, in the future, advances in the measurement of LAI at the local scale (LIDAR) will allow finer model calibration and validation.

1160

1165

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

**Figure 14 | Projected (by PHOREAU) against observed satellite leaf area index (LAI) for all 340 RENECOFOR and ICP II validation inventories.** The y-axis shows the LAI predicted by the model from the stand inventory at the start of the simulation, while the x-axis represents the PROBA-V LAI value for the machining coordinate and inventory year, averaged between July, August and September. Stand points are color coded by dominant species (see legend in bottom left). The size of points shows inventory basal area. The dashed red-line is the 1:1 line; the black full line represent the regression line of the linear model between observed and predicted LAI, with confidence interval represented by the grey shaded area. Associated statistics in Table S6.

1170

1175

## 4.3 Evaluation against tree basal area increment

PHOREAU demonstrated satisfactory predictive capability for tree-level mean annual basal area increment (BAI) across diverse species and climatic conditions throughout Europe (Fig. 15). The model achieved a strong correlation between observed and predicted values ($r$ = 0.68, $p$ < 0.001, n = 81655; Table S4), with satisfactory levels of prediction accuracy (RMSE = 0.00106, AB = 0.225, and AAB = 0.793). However, the model dampened

1180
<figure><img src="image_103.png" /></figure>

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

the observed variability in tree growth, tending to underestimate the productivity of the most vigorous trees while simultaneously overestimating growth of the least productive trees.

When assessed at the species level, the Pearson correlation coefficients varied substantially, from 0.14 for *C. avellana* to 0.913 for *U. glabra* (Table S4). Prediction accuracy also differed widely among species, with an average RMSE of 0.00103 and an AB of 0.34. Correlation coefficients were generally higher for the 13 main species of the study (those that dominate at least one of the 340 simulation inventories) compared to secondary species (average *r* = 0.60 and 0.53, respectively), with a pronounced tendency for the model to underestimate the productivity of these secondary, generally understory species, whose growth rates were not recalibrated on forest growth data in the ForCEEPS study (Morin *et al.*, 2021).

In comparison with the ForCEEPS model, which was applied to the same dataset (Fig. W13, Table S4), PHOREAU demonstrated a moderately improved performance in predicting tree productivity. It yielded higher Pearson correlation coefficients, as well as lower RMSE and absolute errors. Despite these improvements, PHOREAU's predictions exhibited a comparatively greater average bias.

<figure><img src="image_105.png" /><figcaption>**Figure 15** | Projected (by PHOREAU) against observed mean annual tree basal increments (BAI) for all simulated trees over the 340 RENECOFOR and ICP II validation inventories. Tree points are color coded by species (see legend above). The dashed red-line is the 1:1 line; other full lines represent the regression lines of the linear model between observed and predicted tree productivity, with confidence intervals represented by the grey shaded area (in black the overall regression; coloured lines for species-specific regressions). Species-specific regressions are only shown for stand dominant species (in bold in legend) Associated statistics for the global simulation in top left, while species-specific statistics can be found in Table S1.</figcaption></figure>

1185

1190

1195

1200
<figure><img src="image_106.png" /></figure>

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

## 4.4 Evaluation against stand basal area increment

At the stand level, PHOREAU exhibited robust performance in reproducing mean annual BAI across most species and environmental conditions. Overall, there was a strong correlation between observed and predicted values across all 340 simulations (r = 0.62, p < 0.001) with a small margin of error between observations and predictions (RMSE = 0.23, AB = 3.7%, and AAB = 0.34; Fig. 16, Table S5). However, the accuracy varied when species were analyzed individually. While the model generally showed no systematic bias (RMSE = 0.23; AB = -2.2% on average), some species exhibited notable biases and variability, particularly in the most productive plots where the model tended to underestimate productivity. This was especially evident for *P. halepensis* (RMSE = 0.35, AB = -65%) and *P. menziesii* (RMSE = 0.18, AB = -29%), though both had relatively small sample sizes. Even for *Q. petraea* (RMSE = 0.19, AB = -17%), where sample size was not a limitation, a similar bias was observed.

When examining the relationship between prediction errors and various stand characteristics (Fig. W11), no strong systematic biases were identified with respect to site-specific factors such as rainfall, temperature, stand density, or simulation duration. However, the regression analysis revealed a weak but statistically significant positive relationship between errors and site water-holding capacity (SWHC) (slope = 0.0034, r = 0.138, p < 0.05), suggesting a tendency to underestimate productivity on drier soils. Additionally, there was a strongly significant negative relationship between errors and initial stand basal area (slope = -0.0044, r = -0.21, p < 0.001), indicating that the model underestimates productivity in the most productive stands.

In comparison to the ForCEEPS model applied to the same dataset, PHOREAU demonstrated enhanced predictive accuracy across all evaluated metrics. PHOREAU produced a higher Pearson correlation coefficient than ForCEEPS (r = 0.62 vs. r = 0.53 respectively), along with lower RMSE (0.23 vs. 0.316), average bias (AB = 3.7% vs. 7.7%), and average absolute bias (AAB = 0.34 vs. 0.44; see Fig. W14, Table S5). These results highlight PHOREAU's improved capability in predicting stand productivity compared to ForCEEPS.
<figure><img src="image_108.png" /></figure>

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

<figure><img src="image_110.png" /><figcaption>Figure 16 | Projected (by PHOREAU) against observed mean annual stand basal increments (BAI) for all 340 RENECOFOR and ICP II validation inventories. Stand points are color coded by dominant species (see legend above). The dashed red-line is the 1:1 line; other full lines represent the regression lines of the linear model between observed and predicted stand productivity, with confidence intervals represented by the grey shaded area (in black the overall regression; colored lines for species-specific regressions). Associated statistics for the global simulation in top left, while species-specific statistics can be found in Table S2.</figcaption></figure>

1230

1235

## 4.5 Evaluation against potential natural vegetation data

When comparing the distribution of predicted dominant tree species after 2,000-year simulations along the environmental gradient covered by 250 sites across Europe (Fig. 17), the model performed well, with 62% of predictions accurately matching observed community compositions, and 24% partially accurate predictions (outperforming ForCEEPS’ 43% accurate predictions). Yet, PHOREAU’s ability to accurately predict potential natural vegetation (PNV) varied depending on site conditions, with a noticeably larger uncertainty for Mediterranean forest types, humid beech forests, and mixed montane spruce-beech forests. A detailed view of the predicted dominant species (Figure W16) revealed that much of this uncertainty stemmed from PHOREAU’s tendency to overestimate the competitive advantage of *Q. robur* relative to *Q. petraea* and *F. sylvatica* in both hot and mild climates. Despite these discrepancies, the model demonstrated strong predictive performance in extreme environments, accurately predicting species composition at both extremely cold and extremely warm sites.

1240

1245
<figure><img src="image_111.png" /></figure>

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

Figure 17 | (a) Distribution of the 250 tested sites in the PNV diagram of supposed dominating species (built according to mean annual temperature and precipitation sum). PNV dominating species are Pc (P. cembra); Pa (P. abies); Aa (A. alba); Fs (F. sylvatica); Qp (Q. petraea); Qr (Q. robur); Pp (P. pinaster); Ph (P. halepensis); Qi (Q. ilex) Circle colors indicate the agreement between simulated and PNV dominating species after the 2000 years PHOREAU simulations. Green: sites for which the dominating species was accurately predicted. Orange: sites for which the second-ranked (by basal area) species was accurately predicted, but not the first-ranked. Red: sites for which neither the first-ranked nor second-ranked species were accurately predicted.

(b) Geographical repartition of the 250 sites (RENECOFOR and ICP II) used for PNV validation, colored by potential niche composition. Shapes indicate prediction success, as described above.
<figure><img src="image_113.png" /></figure>

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

1250 **5 Discussion**

**5.1 A process-based model to investigate diversity-productivity and diversity-resilience relationships**

The difficulties inherent in integrating trait-based processes in a semi-empirical framework justified evaluating PHOREAU on a variety of metrics — including predicted foliage area, soil water and stem water potentials — which, to our knowledge, has never been attempted before, at least for this kind of model. Furthermore, the bottom-up approach we have adopted mitigated the risk of error compensation and of equifinality, which often appear when some parameters or processes covariate and compensate each other in respect to an integrative metric. Avoiding equifinality was crucial to the development of PHOREAU, because as climatic conditions deviate from the historical baseline in future years, correlations between processes that were equifinal for historical conditions may shift, limiting the ability of the model to accurately predict the impact of climate change on forest functioning. While direct validation on annual growth is rarely done for gap models because of the inherent difficulty of reproducing such metrics for models not originally designed to work at such short temporal scales (Mette et al., 2009; Fyllas et al., 2014), the more granular representation of stand functioning of PHOREAU justified our evaluation on short-term individual tree and stand productivity. The good performance of the model across the wide range of species and conditions used in the productivity and PNV validation — including Mediterranean and boreal forests — demonstrates its widespread applicability to European forest ecosystems. Furthermore, the state-of-the-art validation dataset used in this study will serve as a baseline to assess any further refinements to the model, as additional species traits become available.

1255

1260

1265

1270

In contrast to ecophysiological process-based models than can be parametrized using only physiological and functional traits (Davi *et al.*, 2005b; Maréchaux and Chave, 2017), PHOREAU eschews a direct representation of carbon assimilation and allocation, in favor of a growth-reduction based approach. While this simplification does distort actual tree functioning and ignores the importance of carbon reserves in buffering year-on-year growth (Körner, 2003), it presents a number of advantages when considering the ecological processes that shape species composition. In addition to a significant gain in computing time, it curtails the uncertainty in model predictions that can result from equifinality, by limiting the number of variables directly impacting growth. Furthermore, by calculating tree growth, leaf area, mortality and establishment rates on the basis of well-established observed parameter values, to which process-based reductors are subsequently applied, we were able to maintain realistic stand basal and foliage areas over the length of the simulation. This result is a prerequisite to any temporal exploration of diversity-resilience relationships in drought-stressed forests: only by accurately predicting the evolution of forest foliage and basal area can we then study the effects of species-mixing (Forrester and Pretzsch, 2015) for forests functioning at eco-hydrological equilibrium. This is why our integrative validation on the ICOS sites is an important milestone in the development of hydrology-based forest models: unlike usual hydrological validations (Morales *et al.*, 2005), not only did PHOREAU provide robust predictions of water fluxes for many years over a diverse set of conditions and species, it did so with no *a priori* fixing of stand leaf and basal area, instead calculating the evolution stand structure on the basis of water-stress feedbacks.

1275

1280

1285
<figure><img src="image_115.png" /></figure>

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

1290

## 5.2 Limitations and future avenues of improvement

Despite good correlations and low average bias, PHOREAU predictions consistently underestimated the observed variability across almost all considered metrics, including soil water quantities, stem water potentials, tree productivity, stand productivity, and stand foliage areas. This attenuating effect is in itself not surprising given the necessary simplifications presented by any modelling approach, and results from a number of unavoidable factors: precision of climatic and soil texture data (especially for ICP II sites); utilization of single sets of species parameters disregarding intra-specific genetic and phenotypic trait variability; lack of 3D representation of competition among trees. While climatic, soil, and species traits inputs can easily be refined for more granular simulations at the local and regional level, taking into account site exposition and fertility, the strong hypothesis of the PHOREAU model regarding the horizontal homogeneity of competition for light and water inside a patch will always be an obstacle to capturing the individual dynamics of trees advantaged or disadvantaged by microtopography and spatial allocation of tree crowns and rooting systems. Despite this inherent limitation, the integration in PHOREAU of many previously disregarded or implicit processes, including explicit roots, phenology, process-based tree hydraulics, and microclimate, has allowed it to outperform the ForCEEPS model in better predicting both short-term growth and long-term species composition. Furthermore, the gap between the two models' predictions is likely to become greater under future conditions, where PHOREAU is expected to be more robust as it explicitly represents key processes, such as drought stress and phenology, in a more mechanistic way.

1295

1300

1305

1310

However, by introducing a more granular representation of tree functioning, PHOREAU has induced a mismatch between some of the parameters used in the model and the role they were originally intended and calibrated for. This mismatch, particularly evident for the optimal species growth rate parameter ($g_s$) and for foliage allometry parameters, is responsible for the difficulty in reproducing the growth of extremely productive trees, and the overall underestimation of the productivity of species like *P. halepensis*, *F. excelsior*, or *A. pseudoplatanus* (see Table S4). Because the optimal growth rate in ForCEEPS was calibrated for the main French species based on the top 10th percentile of annual diameter incements measured in the NFI database (IGN, 2020) and for other species dates back to even earlier studies (Didion et al., 2009), it is in reality more akin to a growth rate under relatively unconstrained conditions than an actual optimum. As we updated the model's representation of light and water use constraints to a more process-based approach, we have likely introduced constraints already implicitly present in this aggregated growth rate parameter, essentially penalizing trees twice for the same factor. As we continue to refine the PHOREAU model, a major challenge will therefore be recalibrating this parameter to better reflect actual potential growth unconstrained by competition, despite inherent difficulties in obtaining such data (Pretzsch, 2009).

1315

1320

1325

Similarly, the parameters with which foliage area is derived from tree diameter have not been fully updated to reflect the new importance of foliage area in driving modelled water fluxes. Despite the many changes introduced in the representation of tree crowns and the partial validation on satellite data, the model demonstrated a poor ability to predict measured litter LAI for sites of similar composition and basal area. Furthermore, neither satellite nor litter-derived total LAI measurements can be used to properly evaluate the predicted vertical distribution of leaf area. However, predicted vertical LAI distribution, from which microclimate and individual light-competition

1330
<figure><img src="image_117.png" /></figure>

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

constraints are derived, is key to model ecological processes, and should therefore next be examined and validated against ground or airborne LIDAR and microclimate measurements.

Another obvious area of improvement for the model will be a deeper integration of the plant phenology component with other modelled processes. In this study leaf unfolding, leaf senescence, and probability of fruit maturation were computed yearly for an average individual of each species. This method captured inter-specific differences in phenology and temporal light partitioning, but did not account for intra-specific shifts in phenology caused by stand structure. By integrating model variables like microclimate, light availability, and water stress as inputs for an individual-based phenology calculation, PHOREAU will be able to captured well-established variations in leaf phenology between trees of different sociological status (Augspurger and Bartlett, 2003; Cole and Sheldon, 2017; Gressler et al., 2015; Schieber, 2012), which are responsible for the persistence of shrubs and saplings in mature forests (Gill, Amthor and Bormann, 1998; Vitasse, 2013).

## 5.3 Applications and future research perspectives

### 5.3.1 Establishing baseline available water: retro-engineering PHOREAU to predict rooting depths

One of the main causes for the model's attenuation of variability in stand and tree productivity was the uncertainty regarding the actual quantity of soil water available to the trees. This uncertainty is itself the result of a twofold gap in information: lack of data for the texture of deeper soil horizons, and the extremely simplified framework used to estimate tree rooting depths. By choosing to reduce the wide observed differences in rooting depths across biomes (Canadell et al., 1996; Schenk and Jackson, 2002; Fan et al., 2017) and species (Sperry et al., 2002; Fan et al., 2017) to a simple equation based only on tree size and an aggregate drought index based on past climatic conditions, we intentionally avoided any integration of model results (such as tree foliage area or percentage of embolism) in the calculation of rooting depths, as this would have resulted in an optimization of soil available water on precisely the variables we were trying to validate. Unlike other process-based models validated on stand hydraulic fluxes (Ruffault et al., 2023), the fact that PHOREAU produced robust multi-year predictions without using observations to control for stand leaf areas, rooting depths, or actual available water, confirms its possible applications to making realistic dynamic predictions across a large range of forests where this data is not available.

To overcome difficulties related to the soil water parametrization, an alternative approach could be used. For instance, based on the hydrological equilibrium hypothesis (EHE), which states that, in a given edaphic and climatic environment, trade-offs between vegetation water use and drought stress drive canopy density and forest composition toward an optimal hydric state (Eagleson, 1982; Caylor, Scanlon and Rodriguez-Iturbe, 2009), and following the well-substantiated hypothesis that trees function near the point of catastrophic hydraulic failure with narrow safety margins (Tyree and Sperry, 1988; Choat et al., 2012), a retro-engineering of PHOREAU could be realized where rooting depths are calculated by optimizing tree available water such that, for a given inventory and soil profile (Kirchen et al., 2017), foliage area is maximized (Grier and Running, 1977), and plant minimum water potentials are constrained to values to the point of catastrophic xylem failure. Compared to similar EHE-based statistical (Nemani and Running, 1989) or process-based (Cabon et al., 2018) modelling approaches, this retro-engineering of PHOREAU will natively integrate many inter- and intra-specific niche and competition
<figure><img src="image_119.png" /></figure>

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

processes that are integral to forests' actual water use. It will furthermore be a necessary first step in establishing a historical baseline when using the model to predict the medium-term impact of global change on forest composition and functioning, as available water is a major determinant in predicting drought-induced die-off events (Allen, Macalady, Chenchouni, Bachelet, McDowell, Vennetier, Kitzberger, Rigling, Breshears, E.H. (Ted) Hogg, et al., 2010; Anderegg et al., 2013; McDowell et al., 2013).

1375

### 5.3.2 Unraveling the effects of trait diversity on competition and coexistence

The novel approach presented in this study, integrating plant functional traits in a forest dynamics model, was developed to improve the generality of the calibration for new species, but also to cope with the difficulties encountered by ecologists when testing hypothesized links between trait diversity, species competition and coexistence. While differences in traits governing resource use should, intuitively, translate into niche differences that maintain coexistence through competition reduction, attempts to directly link trait dispersion with historical species coexistence have proven challenging (McGill et al., 2006; Adler et al., 2013). This challenge arises from the fact most traits impact competition for several resources at the same time, and that even a temporary advantage in growth can actually result in a lower global fitness when considering population dynamics, with for example feedbacks on drought-induced mortality (Forrester and Pretzsch, 2015) or frost damage due to early onset leaf unfolding (Bigler and Bugmann, 2018). To overcome this difficulty, process-based models of resource competition with processes explicitly relying on species traits have been proposed as a way to unravel the mechanisms linking trait diversity to forest functioning (Levine et al., 2024). Because the effects of climate change on forests will likewise be mediated by complex species mixing effects, the need to develop mechanistic models that bridge the gap between trait-based and ecology and empirical modelling has become urgent to assess the short and medium-terms effects of global warming on existing forests, and discriminate between the possible management scenarios available to forest managers.

1380

1385

1390

1395

The PHOREAU model, having been directly evaluated for most of its processes, could be used as a relevant tool to identify thresholds conditions for species coexistence, dominance, or extinction. A first parsimonious approach could simply consist in identifying the main processes — phenology, water-use, or competition for light — limiting a species fitness at the edges of its predicted distribution (Morin, Augspurger and Chuine, 2007). A more involved exploratory protocol could follow the methodology outlined in Levine et al. (2024). By considering predicted species compositions for a wide range of climatic and edaphic conditions, and taking care to distinguish, for each set of condition, the different mechanistic processes which make up a species' competitive fitness, we could establish relationships between aggregated model metrics (for example growth reductors) and underlying species traits. These relationships could then be used to predict the impacts of climate change on forest composition. In parallel to this approach, and as a prerequisite, predicted species compositions should be compared to actual observed compositions, albeit for a much greater set of points than those for the potential composition validation presented in this study, dissipating any remaining uncertainties regarding the representation of regeneration and mortality, which is one of the main current challenges for forest modelling (Cailleret et al., 2017; Vanoni et al., 2019).

1400

1405

1410
<figure><img src="image_121.png" /></figure>

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

### 5.3.3 Evaluating management policies under future climate scenarios

A further policy-relevant application of the PHOREAU model in the coming decades lies in its ability to simulate forest management scenarios under different climate trajectories, and evaluate their outcomes based on key ecosystem service metrics, including wood production, biodiversity conservation, and carbon sequestration. As forests play an increasingly critical role in helping countries meet sustainable development goals (Chapin III et al., 2008), and with forests storing roughly half of terrestrial carbon (Friedlingstein et al., 2019), predicting forest carbon dynamics and its response to management decisions under climate change has become an essential consideration for forest managers. However, while policy makers — supported by the recorded increase in the European forest carbon sink in the early 21st century (Pan et al., 2011) — table on a continued increase in the share of carbon emissions removed by forests (with a target of 40% in France by 2050), this dynamic has already shown signs of slowing (McDowell et al., 2020) as the early forcing effect of climate warming on forest productivity is now counterbalanced by increased drought-induced tree mortality (Allen, Macalady, Chenchouni, Bachelet, McDowell, Vennetier, Kitzberger, Rigling, Breshears, E. H. (Ted) Hogg, et al., 2010; Hammond et al., 2022). While previous studies have evaluated the performance of different management strategies for carbon sequestration over the next decades based on a priori global forest biomass trends and management rules (Bastick et al., 2024; du Bus de Warnaffe and Angerand, 2020), very few models, to our knowledge, have attempted the dynamic integration of forest management with stand-specific future conditions to predict the evolution of the forest carbon stock. By integrating management, growth, and hydraulic processes, PHOREAU is uniquely positioned to simulate more realistic and agile forest trajectories, and to help forest managers by giving them insights about how to better adapt forest to new environmental conditions through management actions.

In conclusion, by combining a detailed representation of plant functional traits with the flexibility required for large-scale simulations and species calibration, PHOREAU offers a unique compromise between ecophysiological realism and operational applicability — making it a valuable tool for both ecological research and forest management under climate change.

1415

1420

1425

1430

1435

1440

1445

1450
<figure><img src="image_123.png" /></figure>

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

# 6 Supplementary information

## Appendix A: Decoupling tree height from diameter: light-dependent plasticity

The predictive power of gaps models is tied with their representation of stand structure. Yet most classic gap models, including ForCEEPS, do not simulate a dynamic tree height, instead inferring it from the tree trunk diameter through an allometric relationship. It follows that for a given species, every individual follows the same height-diameter trajectory. While this is consistent with the fact most forestry surveys report basal diameter without height, this means that the models cannot represent site effects on maximum height, as well as the effects of competition for light on the height-diameter relationship. In reality dominated understory trees tend to allocate more carbon to height growth than diameter growth. Conversely, trees in low-density or thinned forests have greater diameter growth and slower height growth (Oliver and Larson, 1996). Furthermore, this sensitivity of growth allocation to competition for light is more marked in shade-intolerant species (Delagrange et al., 2004).

The effects of competition for light on growth allocation are crucial for understanding stand dynamics, as small initial differences in height tend to increase with time unless corrected by greater height growth. Forest managers have long known that tree maximum height varies from site to site with tree age and density (Fortin et al., 2019), and forest growth models often use different height-diameter depending on site conditions (Mehtätalo, Miguel and Gregoire, 2015). Attempts to implement dynamic height growth in gap models have been shown to increase the realism of simulated stand structure, without reducing general applicability. For instance Rasche et al. (2012) have implemented such a dynamic height in the ForClim model on which ForCEEPS is originally inspired. Instead of the static relationship between diameter and height (h), height increments are calculated at each time-step $\Delta H = f_h \Delta D$ through a function $f_h$ that distributes growth between diameter and height growth according to a competition-for-light driven parameter s, which replaces the original fixed species-specific allometric parameter. Since the yearly diameter increment uses previous-year height in its calculation, its formulation also had to be adapted to account for the fact that height is dynamic and no longer directly calculated from diameter. These adaptations have been used in our modified ForCEEPS model, albeit with two important modifications.

Firstly, the parameters of the growth-distribution coefficient $g_s$ were adapted to be more conservative, and better reflect the species-specific relationship that had already been parametrized:

$$f_h = s \times \left(1 - \frac{H-137}{H_{max,s}-137}\right) \qquad \text{Eq. A1}$$

$$s = s_{original} + (kLa * 10) * (1 - AL_H) \qquad \text{Eq. A2}$$

where $kLa$ is the species shade-tolerance, $H$ the tree height in centimeters, $H_{max,s}$ the maximum species height, and $AL_H$ the light availability at the top of the tree crown.
<figure><img src="image_125.png" /></figure>

Secondly, we adapted the yearly growth equation. In the original formulation by Rasche et al. (2012), because yearly growth is calculated on the basis of total diameter at the start of the year, a tree that allocated more growth to height than to diameter due to competition in year n would have less total growth for year n+1 than a tree that had allocated more growth to diameter, all else being equal. This is a result of the simplifications of the ForClim model, in which the diameter increment is calculated on the basis of previous year diameter instead of the previous year volume. This means tree biomass is only dependent on tree diameter, disregarding its height. This effect has major implications, as originally taller but thinner trees end up with smaller final height and diameters than in the original formulation. A possible solution would have been to replace trunk diameter by volume in the growth equations; but this would have meant reshaping the model from the ground up, and making it less applicable to classic forestry datasets, as actual volume data are very rarely available. In the end, we adopted an *ad-hoc* solution by giving each tree two sets of heights and diameters : a static set ($D_{static}$ and $H_{static}$), calculated from the old equations and static allometry relationships, that were only used as an *ad-hoc* proxy for real tree volume in the updated diameter increment equation (Eq. A3) and the calculation of slow-growth mortality (to avoid killing off trees that allocate too much growth to height); and a real set ($D$ and $H$), using the updated equations and dynamic allometry, that was used in all other cases including the light-competition module.

$$ \frac{\Delta D}{\Delta t} = kG * D_{static} * \frac{\left(1 - \left(\frac{H_{static}}{H_{max}}\right)\right)}{2 * H_{static} + f_{h} * D_{static}} \qquad \text{Eq. A3} $$

## Appendix B: Crown-length reversion

The dynamic change of tree crown length was modified to better represent the feedbacks between stand structure and competition for light. In PHOREAU, light availability impacts growth directly and indirectly: directly through the shading growth reduction factor, and indirectly through the crown-length growth reduction factor, which represents long-term crown shrinking due to shading. Individual tree crown lengths are calculated as the product of tree height, and a variable ratio that depends on species characteristics and tree status. This ratio changes according to the light exposition of the tree, between two extreme species-specific values as described in Morin et al. (2021). In the original ForCEEPS framework, seedlings started with a crown ratio set at the species maximum, which then decreased over the tree's lifetime with shading. In particular, this formulation assumes that the crown ratio can only ever decrease or stay the same from one year to the next, with no possibility of reversion when more light becomes available.

Therefore, we have implemented the possibility of crown ratio reversion in PHOREAU. A constantly decreasing crown ratio assumes no increase in light availability over a trees lifetime, disregarding the impact that the death or removal of one tree can have on its neighbours by enhancing light availability and leading to larger crown sizes and denser canopies (see Juchheim, 2020, and Saarinen et al., 2022). We have consequently adapted the original ForCEEPS crown ratio equation to reflect this, with a yearly increase capped at 5% of the difference between the previous-year crown ratio, and the potential crown ratio given current light availability. We are aware this approximation does not take into account the fact that younger trees recover their crowns better due to having more remaining growth potential (Hynynen, 1995).

1495

1500

1505

1510

1515

1520

1525

1530
<figure><img src="image_126.png" /></figure>

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

## Appendix C: *Species-dependent crown shapes*

An accurate representation of crown shapes is an integral component to any model of light competition and canopy interactions between trees (Krůček *et al.*, 2019). In reality the crown shape of any given tree is a complex combination of genetic, allometric, and environmental factors, as crown shape varies across species, age groups, climate, local conditions and the shading status of the tree (Oliver and Larson, 1996). Canopy packing in mixed forests can be partly attributed to this heterogeneity and plasticity of crown shapes, as trees suffer relatively less competition for a given foliage density (Longuetaud *et al.*, 2013).

1535

Crown-shape representation in PHOREAU iterates on the ForCEEPS framework, which already allowed for stratified distributions of foliage area over a vertical axis (Morin *et al.*, 2021). Compared to the previous iteration, PHOREAU allows trees to have other crown shapes than the default inverse-cone – such as conical or ellipsoidal shapes. This is meant to represent broad patrons in crown geometry observed at the European Scale, such as the fact species present in higher latitudes or latitudes tend to have more columnar or conical crowns to capture light coming from a perpendicular angle, whereas species as lower latitudes are more frequently flat-topped for maximum exposure (Kuuluvainen and Pukkala, 1989).

1540

1545

While the lack of explicit tree positions prevent PHOREAU from recreating the asymmetrical crown shapes which result from horizontal constraining between crowns (Niklaus *et al.*, 2017), this simple approach allows for a more accurate representation of side-shading between trees, and captures the way shaded trees tend to become more flat-topped as they reduce their crown height (Oliver and Larson, 1996), while saving some simulation time. See Figure 7 for a visualization of the new crown shapes.

1550

## Appendix D: *Density-dependent light availability*

Any representation of forest canopies and light dispersion has to strike a balance between predictive power — how much photosynthetically active radiation (PAR) does a given tree actually receive at a given moment in time? — and computing cost: by aggregating leaves on a tree-by-tree basis and disregarding differences in angle and light absorption between sun and shade-leaves (Givnish, 1988), by calculating at yearly time-step, and by considering only the vertical stratification without an explicit representation of trunk distribution across space, ForCEEPS is able to compute in a timely fashion what would otherwise take orders of magnitude longer with a more bottom-up approach from the leaf to the tree.

1555

1560

PHOREAU does not diverge from this general framework, which is well suited to working on large-scale inventories (that usually come without tree-level coordinates), and does not suppose any *a priori* knowledge on canopy composition. However, this simplification is not without its drawbacks. Because the light availability of a given canopy layer depends solely on the foliage area present in the layers above it, with no accounting for how this foliage is actually distributed, light competition is — in effect — boiled down to a single value: the LAI. Intuitively we understand that this does not quite tally with reality: two superposed leaves will intercept less light, all else being equal, than two leaves on a level plane; forests are not horizontally homogeneous, and gaps in the canopy may form as trees die off, allowing saplings to sprout and grow even in dense stands (Nicotra, Chazdon

1565

1570
<figure><img src="image_128.png" /></figure>

and Iriarte, 1999). Due to the links between patchy structures of light availability and tree species diversity and coexistence (Moora *et al.*, 2007), measuring and quantifying microsite light availability has been a focus of research (Parent and Messier, 1996; Tymen *et al.*, 2017), with important implications for forest management (Coates *et al.*, 2003).

This structural limitation — which can be important, e.g. to accurately predict species richness in relation to management — can never be fully worked around. And, in keeping with the general philosophy of the model to strike a balance between complexity and genericity, we opted not to incorporate a complex 3D tree-level light absorption model (le Maire *et al.*, 2013). However, in the transition from ForCEEPS to PHOREAU, some steps have been taken to at least partially account for the horizontal stand structure. This was done in an indirect way by using information available to the model: the stand density.

As in most gap-models, foliage area in ForCEEPS is translated into light availability using a modified logarithmic Beer-Lambert law, see Eq. D1, where light availability is a function of foliage area and a light extinction coefficient $\lambda$. In the original formulation of the law this extinction coefficient is calculated by integrating over the path of the light ray the absorbance and density of the materials it crosses. This calculation — which accounts for the angle of the leaves, the angle of the sun's rays, the different absorbances between species and sun and shade-leaves, and the distribution and clumping of the leaves and trees (Smith, 1993; Dufrêne and Bréda, 1995) — is usually simplified into an empirical constant extinction parameter, which can vary from site to site (Vose *et al.*, 1995; Binkley *et al.*, 2013). However, in the ForCEEPS framework, where stand composition is an emergent property and not an input, a single $\lambda$ value is used regardless of site conditions.

Following the methodology outlined in (Nilson, 1971; Black *et al.*, 1991; Bréda, Soudan and Bergonzini), PHOREAU integrates a clumping factor $\Omega$ in its calculation of the light extinction coefficient. This clumping factor ranges from 0 (corresponding to a fully concentrated distribution) to 1 (corresponding to a perfectly homogeneous distribution), and represents the aggregation of leaves within each tree and between the trees themselves. The advantage of this approach is that $\Omega$ can be calculated each year as an emergent variable, allowing the model to capture observed trends like the inverse relation between LAI and the light extinction coefficient (each additional increment of leaf area blocks marginally less light) (Dufrêne and Bréda, 1995). The clumping factor in PHOREAU is calculated using Curtis relative density (Smith, 1993; Curtis, 1982): with this formulation (see Eq. D2) for a given LAI, a dense stand with small trees will block out more light than a stand populated by a few large trees. This approach is similar to the one used in LAI estimation with MODIS or hemispherical photography, where clumping indices are also used to correct the raw measured LAI (Demarez *et al.*, 2008; Chen *et al.*, 2012; Zhu *et al.*, 2018).

A further step would be to incorporate species-specific absorbance values, as leaves of different species react differently to incoming light (Binkley *et al.*, 2013), but this would necessitate gathering data at the species level (data which is, to our knowledge, available only for a select few species). Another possible refinement would be to incorporate the angle of incoming light in the calculation of light availability (Smith, 1980); but this would require modifying the light competition calculation to consider site effects related to slope and exposition.
<figure><img src="image_129.png" /></figure>

$$Light Availability_i = e^{-\lambda * \sum_{j=i+1}^{n} (FoliageArea_j)} \qquad Eq. D1$$

$$\lambda = 0.25 * \frac{BasalArea}{\sum_{j=1}^{N} (DBH_j^2)} \qquad Eq. D2$$

n: total number of layers ; i: layer rank ; N: number of trees

## Appendix E: Incorporation of Specific Leaf Area

The relation between trunk diameter, crown biomass, and foliage area in ForCEEPS are governed by a set of simple allometric relationships calibrated for a few of the main temperate European species, using experimental data collected in Switzerland by destructive sampling in the 1940s and 50s (Burger, 1951; Bugmann, 1996). The refinements that ForCEEPS implemented regarding crown plasticity and explicit vertical stratification were built upon this foundation but did not challenge its underlying assumptions (Morin et al., 2021). This became problematic as the model — and PHOREAU in particular — incorporated more species from a larger geographic range: understory or Mediterranean species in particular that were not represented in the initial calibration dataset. This was directly reflected in model predictions, for example with an overestimation of Quercus ilex or Pinus halepensis mortality due to inflated foliage areas.

A simple solution to this issue was implemented by recalculating the c₂ parameter (used in ForCEEPS to derive a tree's foliage area from its diameter) using a specific leaf area (SLA) value for each species. The retained SLA — the surface area for a given mass of leaves — are those of average adult individuals of each species over a large set of sites (Kattge et al., 2020; Devresse et al., 2024). This new formulation (see Eq. 17) allows the model to capture inter-specific differences in drought resistance strategies (Greenwood et al., 2017), while disregarding for the moment SLA plasticity to tree age, competition, and site conditions (Gratani, 2014).

## Appendix F: Microclimate derived from stand-structure

By integrating fine hydraulic and phenological mechanisms in the overall framework of a forest-structure gap model, PHOREAU has the opportunity to capture the effects of microclimate on plant functioning. Because forest canopies absorb or reflect the majority of incoming solar radiation, reduce wind speeds, convert solar energy into latent heat through evapotranspiration, and block outgoing infrared radiation, climatic conditions in the understory are often buffered compared to those at the top of the canopy, with cooler more stable temperatures during the day, and warmer temperatures during cold nights. This climate dampening effect is more marked for temperature extremes, and for tall, structurally complex dense canopies (De Frenne et al., 2021). Furthermore, it is an important factor in ability of young, understory trees to resist droughts despite their shallow root systems (Forrester and Bauhus, 2016). Because PHOREAU evaluates drought-stress at an individual level by calculating tree fluxes, it can easily make use of microclimatic data for temperature, air humidity, and light availability, to better compute plant evapotranspiration and in turn differentiate water stress among individuals of different heights. In addition, because PHOREAU simulates many small patches each sharing a soil and a canopy height profile, the incorporation of microclimate could help the model capture forest landscape mosaic dynamics, where forests with
<figure><img src="image_130.png" /></figure>

heterogeneous patches are able to host more diversity due to differentiated microclimatic effects on regeneration and drought (Pincebourde et al., 2016).

To derive microclimate temperature and air humidity from macroclimate, we implemented a version of the statistical model developed, calibrated and validated in Gril et al., (2023) and Gril, Laslier, et al., (2023). This model, which has the advantage of using only easily available patch characteristics, uses a simple slope and equilibrium approach, presented in Figure S1, to compute microclimate temperature at soil level ($T_0$) from the corresponding hourly or daily macroclimate temperature ($T^j$). The slope ($m_{slope}$) captures the linear relationship between microclimate and macroclimate, while the equilibrium is the point at which microclimate is equal to macroclimate (Eq. F3). In our case, month mean temperature ($T^m$) is used as the equilibrium. The slope, which acts as a buffer if is lower than 1, is computed daily using patch-level leaf area index (LAI), maximum tree height ($h_{max}$), and vertical complexity index (VCI), as seen in Eq. F4 with corresponding coefficients calibrated over a large dataset of microclimate measurements (Gril, Laslier, et al., 2023). VCI is obtained following Van Ewijk, Treitz and Scott, (2011) by calculating the weighted logarithmic average of foliage area proportion per patch canopy layer ($p_i$), normalized by the total number of layers $n$, as shown in Eq. F5 and Eq. F6. Finally, for any given tree height $h$, the corresponding microclimate temperature $T_h^j$ is derived from soil microclimate and macroclimate using a linear interpolation, as shown in Eq. F1 and Eq. F2.

$$T_h^j = T_0^j + (1 - (w(h)) \times (T^j - T_0^j)) \qquad \text{Eq. F1}$$

$$w(h) = \frac{(h_{max} - h)}{h_{max}} \qquad \text{Eq. F2}$$

$$T_0^j = T^j \times m_{slope} + T^m \times (1 - m_{slope}) \qquad \text{Eq. F3}$$

$$m_{slope} = e^{(0.39 - 0.04 \cdot LAI^j - 0.2 \cdot VCI^j - 0.07 \cdot h_{max})} \qquad \text{Eq. F4}$$

$$VCI^j = - \frac{\sum_{i=1}^{n} p_i \ln(p_i)}{\ln(n)} \qquad \text{Eq. F5}$$

$$p_i = \frac{\text{FoliageArea}_{\text{layer } i}}{\sum_{i=1}^{n} \text{FoliageArea}_{\text{layer } i}} \qquad \text{Eq. F6}$$

j: day or hour; m: month; i: canopy layer; n: number of canopy layers

Calculated hourly microclimate temperatures are then used to compute corrected local vapor pressure deficits (VPD) used in PHOREAU transpiration computations. These temperatures are also used in GDD calculations (see Eq. 26), as well as for seedling establishment constraints based on minimal temperatures ($W_{Tmin}$). For seedlings, soil-level microclimate temperature is directly used; for established trees, the microclimate temperature is calculated the weighted average height of their foliage area distribution.
<figure><img src="image_131.png" /></figure>

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

Because leaf unfolding and senescence dates are integrated in the calculations of LAI and VCI, the slope of microclimate buffering or amplification can change throughout the year.

While this approach presents a number of advantages, it comes with major simplifications. The most important one is certainly the linear interpolation of microclimate over the height of the stand, which neglects actual wind movement and radiation attenuation dynamics. Microclimatic data, measured at different heights below the canopy, would be needed to calibrate a more realistic non-linear function. Other simplifications include disregarding the effect of soil moisture, ignoring horizontal heterogeneity within patches, and assuming monthly mean temperatures are a good indicator of equilibrium.

<figure><img src="image_133.png" /><figcaption>**Figure F1** | Schematic representation of the slope and equilibrium microclimate approach, reprinted from Gril, Laslier, et al., (2023).</figcaption></figure>

**Appendix G: Treewise aggregation**

Because the runtime of a SurEau simulation is driven by the number of distinct water-holding compartments — the atmosphere, soil layers, and mostly importantly tree organs — the first step to reducing the runtime of a SurEau simulation is to reduce the number of initial trees. This approach requires that the global stem volumes and foliage areas remain the same at the stand level, as these are the main drivers of water-use in SurEau and *in natura* (Wullschleger, Meinzer and Vertessy, 1998). The aggregation method ensures this through by summing and averaging, at the cost of some precision in the description of the competition for water.

The degree of simplification is specified at the start of the PHOREAU simulation by choosing a number of *classes*: this is the maximum number of aggregate trees created per species at the start for each SurEau run-year. It follows that, for example, a three-class aggregation in a stand with 4 species will result in SurEau initializing with at most 12 trees, which is a more manageable number. To preserve the overall structure of the stand, trees are distributed within classes on the basis of trunk diameter: for an *n*-class aggregation, for each species, the range of diameters between 7.5 cm and the largest diameter at breast height is decomposed between *n* − 1 segments of same size: classes are then created by grouping all the trees with a diameter at breast height located between the extremities

1690

1695

1700

1705

1710
<figure><img src="image_134.png" /></figure>

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

1715 of a given segment, and the last class is composed of all the juvenile trees smaller than 7.5 cm. A consequence of this method is that a class may contain no tree for a given year, and that trees may move between classes from one year to the next as they grow in size.

1720 After the distribution, a single aggregated tree is created for each class. The volume of this aggregate tree is the sum of the volumes of all the trees in the class; its height the average of their heights; its foliage area the sum of their foliage areas; its root depth the average of their root depths; its root biomass the sum of their root biomasses; and finally its light availability the average of their light availabilities. See Figure 5 for an example case.

## Appendix H: Dry-year selection

1725 The second optional way of optimizing PHOREAU performance revolves around modifying the rate at which SurEau is called from ForCEEPS. By default, the two submodels are run on a 1-to-1 basis, with SurEau being called at the beginning of each year; but a more parsimonious approach is to run SurEau only for the driest years of the simulation. This simplification is based on the idea that the impact of drought on forested stands, and especially on tree mortality, does not follow a linear curve, but rather depends on climate extremes, physiological thresholds and tipping points (Hartmann et al., 2018). Because this approach requires a prerequisite ranking of all of the years of the simulation according to their dryness, we use an integrative Drought Index calculated for each year (Morin et al., 2021). The rate of SurEau calls — every two years, five years, etc., — is set by the user before the start of the simulation, with a trade-off between runtime and the accuracy of drought-response predictions. At the start of the simulation, the driest year among the first *n* years is selected as the year SurEau will be called; then, at the start of the *n* + 1 year, the driest year among the next *n* years is selected, and so on.

1740 In assessing the effects of drought events on trees, PHOREAU distinguishes between short-term adaptations and long-term non-reversible consequences — respectively feedbacks on growth and on mortality. The independence of these two mechanisms is key to avoiding confusion between two sources of mortality: that caused by long-term carbon starvation — represented in PHOREAU by diameter growth falling under a certain threshold — and that caused directly by extreme drought through high level of xylem embolism leading to hydraulic failure (Cochard et al., 2021b). A tree subjected to consecutive years of water stress may maintain its conductive vessels but die off due to a lack of carbon intake and defoliation; another may die following a single month of acute water stress despite strong carbon reserves. By establishing a clear distinction between these two pathways, PHOREAU is able to account for the different drought response strategies observed among species.

1745 In PHOREAU, the impact of drought on growth is assessed using the degree of stomatal closure, converted into a drought index *Dr1*. Compared to the original ForCEEPS formulation which uses a simple monthly water budget (Bugmann and Solomon, 2000), this new mechanism takes advantage of the detailed hydraulic framework of SurEAU to account for competition for water as well as inter-specific differences in dealing with water-stress. For seedling establishment — for which SurEAU cannot be used — the original drought index *Dr1* remains used as a proxy for global stand water availability.
<figure><img src="image_136.png" /></figure>

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

Schematically, as soil water reserves become depleted and soil water potential decreases, trees adapt their conductance by closing off stomata in order to reduce water loss and maintain twig and leaf potentials above cavitation thresholds (Cochard, Bréda and Granier, 1996; Cochard et al., 2002). This regulation mechanism prevents the premature death of branches and trees due to uncontrolled embolisms, as trees reduce their water loss until only cuticular transpiration remains. The relation between leaf water potential and stomatal closure is an important trait describing a species' response to drought: constrained by a trade-off between carbon gain and risk of hydraulic failure (Brodribb et al., 2003; Venturas et al., 2018), it is correlated with the more often measured turgor loss point (TLP) (Brodribb and Holbrook, 2003). While the link between turgor loss and reduced growth is well-documented (Cabon et al., 2019; Peters et al., 2020; Potkay et al., 2022), for PHOREAU stomatal aperture was selected as a continuous variable allowing for a finer feedback.

Stomatal aperture $\gamma$ in PHOREAU is derived at each time-step from leaf water potential $P_{L,sym}$ using a sigmoid curve described by two species-specific traits: $P_{gs12}$ the water potential causing 12% stomatal closure, and $P_{gs88}$ the water potential causing 88% stomatal closure (Cochard et al., 2021b). Actual stomatal conductance is then calculated as the product between this stomatal aperture ratio and a maximal stomatal conductance value for a given climate. To calculate the drought reduction index $DrI$ of a given tree, daily stomatal apertures ratios $\gamma_j$ are averaged over the photosynthetic period, which are then averaged over the year (Eq. 13).

## Appendix J: Drought feedback on defoliation

Between the normal closing and opening of stomata to regulate water flow, and the runaway embolisms responsible for tree mortality after prolonged extreme droughts, trees exhibit a range of intermediate responses to water stress. Among these regulatory mechanisms, the adaptation of leaf area to moderate water stress is of particular importance for any model, such as PHOREAU, which integrates tree growth and drought-resistance.

Water limitation impacts leaf area through three main pathways: the premature shedding of leaves, the disruption of new bud formation (Bréda et al., 2006), and plastic biomass allocation to leaves (Martínez-Vilalta, Sala and Piñol, 2004). These mechanisms function at gradually longer time-frames: a cohort of trees may shed their leaves one year in response to extreme drought, and recover their full canopy the next; another may experience several years of decreased leaf area while its leaf phenology cycle is disturbed; and yet another cohort may have permanently shifted to produce less leaf area to adapt to chronic soil water limitations (Limousin et al., 2012; Martin-StPaul et al., 2013). This graduated temporal response is complicated by the fact it is differentially applied among species, following the classic split between drought-avoidance and drought-resistance strategies: indeed, there is evidence that while the reduction of leaf area improves resistance to moderate drought events, it may not avail against severe water stress (Limousin et al., 2022). Furthermore, the short-term gain in drought-resistance of a reduced photosynthetic surface may eventually offset by the negatives consequences of reduced carbon uptake (Poyatos et al., 2013), and the link between leaf area and a reduction of fine root biomass (Gieger and Thomas, 2002).
<figure><img src="image_138.png" /></figure>

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

While the integration of defoliation has been shown to improve the predictions of tree mortality models (Dobbertin and Brang, 2001), this integration is complicated by the fact that few are able to account for the dual role of leaves in carbon-assimilation and water-use. However, unlike most mortality models, the PHOREAU model has the major advantage of being able to disentangle the contradictory effects of leaf area on growth and drought resistance, and of having an explicit representation of the root compartment with water uptake driven by fine roots and ultimately leaf area (see Sect. 3.2.7).

**Appendix K: Drought feedback on mortality**

Drought-induced mortality in PHOREAU is derived from the percentage of cavitation, i.e. the percentage of loss of conductance (PLC). This mortality mechanism is entirely distinct from the pre-existing slow-growth mortality in ForCEEPS, and the previously described drought feedback on growth. Indeed, contrary to the slow-growth mortality that reflects carbon starvation and the long-term integrative effects of dehydration coupled with temperatures and competition for light on the capacity of trees to grow and survive (Bugmann and Solomon, 2000), this feedback is only intended to capture catastrophic water failure caused by extreme drought events, irrespective of the overall prior health of the tree. Unlike the stomatal closure used in drought feedback on growth, the cavitation of a tree's hydraulic system is neither quickly reversible, nor does it follow a linear response to hydraulic stress. Furthermore, it occurs only after the stomata have been closed, when, under extreme stress conditions, residual water flow through the cuticle empties the plant's water reservoirs. As water is drained from the soil and the water potential of the system becomes more and more negative, the conductance of a tree's hydraulic system may remain stable until a certain point is reached, when it rapidly decreases as the xylem vessels are embolized and air are formed (Tyree and Sperry, 1989). This non-linear, tipping point response of conductance loss to decreasing water potentials is described by the *vulnerability curve* of the species. This curve, in the shape of an inverse sigmoid function, is described for each species using a $P_{50}$ parameter. This parameter, responsible for the main differences in drought-resistance between species (Delzon and Cochard, 2014), is the water potential causing 50% cavitation in the xylem (Cochard *et al.*, 2021b).

**Appendix L: The rain interception module**

Capitalizing on the capacity of PHOREAU to predict individual-tree daily foliage area values that integrate allometry, competition, frost, phenology, and drought-defoliation effects, we implement a rain interception module that reduces incoming rain based on the daily leaf area of the stand. Modelling rainfall interception — defined as free water that evaporates from the leaves and barks of trees after a rain event — is an important component for any model trying to water cycles and tree water balance (Granier *et al.*, 1999; Davi *et al.*, 2005a). The intensity of the interception has been shown to grow linearly with leaf area, for values ranging from 20% to 35% of cumulated rainfall in temperate and continental climates (Bréda *et al.*, 2006). While secondary factors such as irradiance, windspeed, and vapor pressure deficit impact the rate of interception *in natura*, as a first approach we have chosen a simple implementation, inspired from Medfate (De Cáceres *et al.*, 2023b), based solely on daily leaf area, rain volume, and potential evapotranspiration.

A canopy storage volume is derived from the foliage area of the stand. This volume is incremented at a daily time-step with incoming rainfall, and outgoing evaporated water. For a given volume of incoming rainfall, the
<figure><img src="image_140.png" /></figure>

throughfall, or the volume of water to reach the ground, is calculated with a simplified Beer-Lambert formula, in a similar fashion to the way light extinction is computed. Because the canopy storage volume is itself limited, any intercepted water that overflows this maximal quantity flows down the soil; a natural consequence of this property is that a given volume of given rainfall will yield a greater cumulated throughfall when concentrated in a single day, than when distributed over several days with intervening evaporation. The algorithm, presented below in Eq. L1, computes the daily stand-wide throughfall volumes that then serve as inputs to the water balance model.

1835

1840

$$ 
\begin{array}{l}
\text{CanopyStorageCapacity}_j = \frac{LAI_j}{2} \\
\text{PotInterceptedRainfall}_j = \text{Rainfall}_j * (1 - e^{-0.5 * LAI_j}) \qquad \text{Eq. L1}
\end{array}
$$

$$ \text{PotThroughfall}_j = \text{Rainfall}_j * (e^{-0.5 * LAI_j}) $$

$$ \text{CanopyStorage}_j = \text{CanopyStorageCapacity}_j - \text{CanopyStock}_{j-1} $$

1845

$$ \text{Throughfall}_j = \begin{cases} \text{PotThroughfall}_j & \text{PotInterceptedRainfall}_j \le \text{AvCanopyStorage}_j \\ \text{Rainfall}_j - \text{AvCanopyStorage}_j & \text{PotInterceptedRainfall}_j > \text{AvCanopyStorage}_j \end{cases} $$

$$ \text{StoredWater}_j = \text{Rainfall}_j - \text{Throughfall}_j $$

$$ \text{CanopyStock}_j = \text{Max}(0, \text{CanopyStock}_{j-1} + \text{StoredWater}_j - \text{PET}_j) $$

1850

j : day of year

## Appendix M: The bootstrap algorithm

In the PHOREAU framework, the leaf area is updated at the end of the year, after each tree's crown length has been updated according to the light availability. However, the light availability that is used to calculate the new crown lengths is the result of the stand area of the previous year, which is itself the result of the previous year's crown lengths. This asynchronicity means that – disregarding other processes like growth regeneration and mortality – the estimation of stand area will oscillate around an equilibrium state. While this equilibrium state is dynamically stable, the oscillations for the first few years are large enough to be significant. This is especially problematic when starting the model from an inventory: because actual crown lengths are rarely available, the model is forced to initiate the crown at the maximum species' value; the resulting very low light availability means that the following year the crown lengths will be reduced by a large factor, which means that more light will be available the year after that, causing a new spike in stand leaf area. It is to correct for this effect that we implemented a bootstrap algorithm where, before the first year of the simulation, multiple iterations of the light competition module are run until the shift in stand area between two successive iterations becomes negligible

1855

1860

1865
<figure><img src="image_141.png" /></figure>

<figure><img src="image_142.png" /><figcaption>**Figure M1 | Illustration of PHOREAU canopy bootstrap algorithm.** Top: one-sided leaf area indices predicted by the PHOREAU bootstrap algorithm, initialized with a Picea abies dominated inventory (RENECOFOR EPC 39a, 2003). Bottom: three snapshots of predicted foliage area and light availability vertical stratification at different steps in the algorithm. For details on the calculation of the Vertical Complexity Index (VCI), refer to Appendix F.</figcaption></figure>

1870

1875

1880
<figure><img src="image_143.png" /></figure>

## Appendix N: The Integrated Carbon Observation System

1885

The Integrated Carbon Observation System is a network of stations that measure ecosystem-atmosphere exchanges of greenhouse gases and energy at a high frequency (Baldocchi, 2003), using the eddy-covariance technique. In addition, a large set of ancillary variables needed for the interpretation of the flux data are also measured: for forest stations these include, among other measures, tree inventories, leaf area index, and soil data — all of which can be leveraged for our modelling purposes (Gielen et al., 2018). The large scope of measured variables in ICOS framework makes any validation based on it easily scalable, and will in the future allow testing of any newly integrated PHOREAU processes (such as carbon retention or vertical micro-climate interpolation). Finally, a set of rigorous specifications for the installation of the eddy-covariance tower sensors, and a common pipeline for the post-processing of the raw data through the Ecosystem Thematic Centre (ETC), ensure the high level of comparability between sites that is necessary for large-scale model evaluation.

1890

1895

## Appendix O: Puéchabon

The Puéchabon experimental site (43°44′30″N, 3°35′40″E, altitude 270 m) is located in a forest of holm oak located in the South of France near Montpellier. With its last clear cut in 1942, and managed as a coppice for centuries before that, the site is characterized by a high density (5000-7000 trees/ha) of small (5.5 meter high overstorey) *Quercus ilex* trees: they make up an old forest with a basal area of 30 m²/ha, (Rambal *et al.*, 2014), and an LAI around 2.2 with little seasonal variability. Located on a flat area, with a rocky soil of Jurassic limestone filled with clay, its small water reserve (roughly 130 mm of water over the 5 meter profile) and typically Mediterranean precipitation pattern (highly variable from year-to-year, with a measured range of 550 to 1550 mm primarily concentrated between September and April) made it an ideal candidate to study the long-term effects of drought.

1900

1905

Within the framework of the *Mediterranean Terrestrial Ecosystems and Increasing Drought* (MIND) project, the diameter of trees contained in twelve 100m² plots have been measured on a year-to-year basis since 2003: these are distributed between three control plots, three thinned plots (33% reduction of basal area), three plots with partial rainfall exclusion (33% throughfall), and three thinned and rainfall excluded plots (Gavinet, J.-M. Ourcival and Limousin, 2019). We have used these plots to run simulations from 2003 to 2020, and assess how the PHOREAU model simulates the effects of tree density on drought resistance.

1910

1915

## Appendix P: Font Blanche

The Font Blanche experimental site (5°40′45″E, 43°14′27″N, altitude: 420 m) is located in a mixed-forest of Aleppo pine and holm oak, with an overstorey of *Pinus halepensis* (13.5 m height) that dominates a coppice of *Quercus ilex* (6.5 m height). With a basal area of 21.3 m²/ha and LAI ranging between 2.5 and 2.7 it is less dense than Puéchabon, but otherwise boasts a broadly similar soil and meteorological profile (Simioni, Marie and Huc, 2016). For our validation we used the 625m² control plot (PM30) of the rainfall exclusion experiment, in addition to the main plot of 6400m² that we split between 25 smaller splots of 267 m² apiece to satisfy PHOREAU homogenous competition assumptions. Our timeframe for this site ranges from 2007 to 2020.

1920

1925
<figure><img src="image_144.png" /></figure>

## Appendix Q: Hesse

The Hesse Experimental site (7°3′59″E, 48°40′30″N, altitude: 300 m) is located in a beech (*Fagus sylvatica*) stand in north-eastern France, on a plain at the feet of the Vosges mountains. Average tree height was measured at 16.2 m in 2005, with a maximum leaf area index over 7.5 indicating a very high level of canopy closure. In comparison to the two previous sites it is characterized by a wetter, semi-continental temperate climate, with a deep loam-clay soil (Davi *et al.*, 2005a; Dufrene *et al.*, 2005). Unlike most sites in the ICOS network it is fertile, fast-growing and subjected to frequent thinnings, with an average tree age of only roughly 40 years in 2005, allowing us to test the capability of PHOREAU to simulate canopy and basal area regrowth after a cut. Furthermore, despite the stand having high rainfall and soil high water holding capacity, droughts events are responsible for most of the interannual variability in tree growth (Granier *et al.*, 2008). We extracted from the inventory four evenly sized 300 m² plots. Because the validation timeframe ranges from 1999 to 2010 when the most data was available (Cuntz *et al.*, 2023e, 2023d, 2023c, 2023b, 2023a) ; ; Betsch et al., 2011; Peiffer et al., 2014; Tuzet et al., 2017; Zapater, 2018), our model also replicates two thinnings that occurred in 2004 and 2009, respectively for 25 and 15% of the basal area.

1930

1935

1940

## Appendix R: Barbeau

The Barbeau experimental site (2°46′E, 48°28′N, altitude: 100 m) is located in the Fontainebleau national forest southwest of Paris. The stand is dominated by sessile oak (*Quercus petraea*) trees that 25 m at 100 years of age, with an understory of hornbeam (*Carpinus betulus*). Mean annual cumulated precipitations of 677 mm are evenly distributed over the year, and feed into a deep soil with roots able to reach at least 150 cm in depth. We initialized our validation over 9 plots of 1000 m² using an exhaustive inventory made in the winter of 2006-2007; we ran running it until 2021, including a thinning in 2011 (Delpierre *et al.*, 2016; Maysonnave *et al.*, 2022). Unlike the other studied sites, growth data was not available on a tree by tree basis, but instead aggregated at the stand level (Briere *et al.*, 2021).

1945

1950

## Appendix S: Supplementary Tables

Tables S1 to S17 are available in the supplements published alongside this article.

1955

## Appendix T: Climate Reconstruction

The SILVAE web portal (Bertrand *et al.*, 2011 and Richard, 2011) offers monthly average temperature and precipitation sum data over France at a finer spatial resolution, accounting for microclimatic differences caused by differences in altitude, exposition, and wind orientation. These time-series, available for the period between 2000 and 2014, were used to correct the coarser ERA-5 Land dataset for all variables except wind-speed: either by direct mean-adjustment for the average temperature and precipitation variables, or after a prior linear regression of the variable over the mean temperature for the given month of the ERA-5 Land time-series. For the average temperature variable, between 2000 and 2014, daily values were corrected by the difference between the average of all the daily ERA-5 Land values for that month and the single monthly value of the SILVAE correction dataset; whereas for the years outside of this range where the corresponding monthly value was not available, the difference

1960

1965
<figure><img src="image_145.png" /></figure>

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

was calculated using the mean of all values for the given month between 2000 and 2014. A similar method was used for the precipitation variable, where the daily values were multiplied by the ratio between summed monthly ERA-5 precipitations, and the single monthly SILVAE value. For the other variables except wind speed the same method was used as for the average daily temperature, except the addition factor was itself first multiplied by the slope of the regression between the temperature and the variable. The wind speed variable was not corrected due to its weak correlation to mean temperature. The workflow for climate reconstruction is summarized in Figure 8.

## Appendix U: Evaluation against leaf area

The importance given to competition for light and leaf area prediction is one of the core principles of ForCEEPS — and the FORCLIM and FORECE models before it. However, because the initial models were focused on long-term forest dynamics, the methodology used to calibrate and validate the light competition module was based on a broad adequation between expected LAI values, and those reconstructed by the model after runs of hundreds or thousands of years starting from the bare ground (Kienast, 1987). Even then, LAI was not usually considered in the final validation, which was made on predicted biomass, basal area, tree density, or species composition (Bugmann, 1996; Wehrli et al., 2006). Notwithstanding the fact that this approach disregards past human interventions in the observed stands, it only accounts for equilibrium states, which becomes problematic when one wishes to apply the model at shorter timescales and consider the shorter-term effects of climate-change on existing forests. Yet, while ForCEEPS did use actual inventories and short-term productivity for its original evaluation (Morin et al., 2021), its performance was not assessed by comparing simulated and observed predicted leaf area index values.

This approach holds up as long as leaf area can be considered to be an intermediary variable. Because the previous models only used leaf area within the framework of their light competition modules, a given tree's predicted leaf area only mattered insomuch as it provided shadow to neighboring smaller trees, decreasing their light availability factor. In this respect, absolute leaf area mattered less than the relative distribution between trees and species, which governed growth and final predicted compositions.

However, in PHOREAU, tree leaf area is also an integral input of another part of the model: the simulation of hydraulic processes. This is because the upwards flow of water through the tree is ultimately driven by the transpiration in the leaves (Ruffault et al., 2022). And, in this respect, water flow is driven not by the relative, but by the absolute quantity of leaf area. Mechanically, a stand with a greater total leaf area index will tend to exhaust its water reserves faster; and tree leaf area, in ecosystems subjected to drought, is directly modulated by recent drought events (Bréda et al., 2006). These mechanisms, which are implemented in PHOREAU, require an accurate prediction of yearly stand leaf area index as a prerequisite condition to any simulation of hydraulic stress.

Unlike other validations of SurEau (Ruffault et al., 2023), the PHOREAU framework prevents the direct use of leaf area index as an input to the model; instead, the model initializes the stand LAI using solely the diameter and height information contained in the initial inventory. This makes the model suited to work on a majority of sites, where trunk diameters are measured but not leaf area, and allows it to make predictions in the future, as the LAI is recalculated on a year-to-year basis. The drawback of this approach is the addition of a new source of error when

1970

1975

1980

1985

1990

1995

2000

2005
<figure><img src="image_147.png" /></figure>

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

2010 LAI is wrongly estimated. This is why, before validating the model on growth or drought-induced mortality, a preliminary validation of the leaf area index predictions was necessary.
<figure><img src="image_149.png" /></figure>

**Appendix V: Height-Diameter Interpolation**

*Height-diameter ratio interpolation.* In order to leverage PHOREAU's ability to reproduce stand light availability

and microclimatic conditions based on the
structure of modelled trees, we used the
newly independent tree height variable (see
Sect. 3.1.2) as an input parameter. However,
height measurements were only available for
a subset of trees across all RENECOFR and
ICP II plots. Therefore, for trees where only
circumference was measured, we applied
plot-specific LOESS local regressions
(Cleveland and Loader, 1996) to estimate
species height-to-diameter curves from
available measurements. The variability in
height-to-diameter relationships among
plots can be seen in Fig. U1 and Fig. S20,

contrasted with the fixed height-to-diameter formula used in the original ForCEEPS framework. The associated
statistics presented in Table S3 highlight the general tendency of the formula to underestimate tree heights in our
study sites (AB = -15.7%; Table S3); this is not necessarily surprising, as the RENECOFOR and ICP II sites
mostly support denser, more productive stands, where trees prioritize height growth to compete for sunlight.

<figure><img src="image_150.png" /><figcaption>Figure U1 | Diversity of site height-to-diameter curves for *Fagus sylvatica*. Refer to Table S3 for details.</figcaption></figure>

2015

2020

2025

2030
<figure><img src="image_151.png" /></figure>

<figure><img src="image_152.png" /><figcaption>2035 Appendix W: Supplementary evaluation figures</figcaption></figure>

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

**Figure W1 | Predicted distribution of stand leaf area and light availability.** This figure illustrates the vertical gradient of predicted light availability indices of the four considered ICOS sites for specific simulation years. The light availability is presented over the aboveground profile, divided into 0.1 m layers. In addition, the area of each shape in the layers represents the predicted aggregate leaf area. Refer to Fig. 5 for light availability index gradient. The figure also includes global annual stand parameters LAI and VCI (see Appendix F of VCI).
<figure><img src="image_154.png" /></figure>

2040

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

**Figure W2 | Predicted fine root area distribution over the soil profile.** For the four ICOS validation sites, for certain simulation years a partial vertical soil profile is shown, with the overall dryness of each soil layer depicted as a gradient using its 10th quantile relative extractable water (REW) percentage. For each species and size class aggregate tree (refer to Appendix G for details on the aggregation method), the distribution of the inverse cone along the soil layers represents the predicted location of its fine roots, with its total aggregate fine root area index (FRAI) shown under. Refer to Fig. W3 for species and cohort color codes.
<figure><img src="image_156.png" /></figure>

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

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

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

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

**Figure W3 | Predicted versus observed evolution of annual stand basal area.** For each simulation site, the bars depict the annual basal area projections generated by the PHOREAU model, broken down by species and size class contributions (refer to Table S15 for associated statistics). The dashed line represents the observed annual total basal area derived from inventory data. Basal area is defined as the cross-sectional area at breast height of all trees per hectare.

<figure><img src="image_161.png" /></figure>
<figure><img src="image_162.png" /></figure>

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

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

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

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

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

**Figure W4 | Predicted versus observed annual stand leaf area index (LAI).** For each simulation site, the bars depict the annual leaf area index projections generated by the PHOREAU model, broken down by species and size class contributions (refer to Table S16 for associated statistics). The dashed line represents the observed annual stand leaf area index (data sources are detailed in Table 1). Leaf area index is defined as the total one-sided leaf area per unit of ground area. Refer to Fig. W3 for species and cohort color codes.

2045
<figure><img src="image_168.png" /></figure>

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

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

**Figure W5** | Predicted versus observed evolution annual stand basal area loss due to mortality. For each simulation site, the bars depict the summed annual total basal area (m²/ha) of all dead trees, broken down by species and size class). Observed values are derived from stand inventories, while predicted values are generated by the PHOREAU model. Also shown are the yearly basal area loss rates, calculated relative to the initial basal area for two distinct time periods in each simulation, along with the total basal area dieback per hectare (G<sub>dead</sub>). Transparent bars indicate years with thinnings (see Appendices Q and R for details), which are excluded from the mortality statistics. Refer to Fig. W3 for species and cohort color codes.
<figure><img src="image_171.png" /></figure>

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

<figure><img src="image_173.png" /><figcaption>**Figure W6 | Predicted versus observed daily real evapotranspiration (ETR).** For each simulation site, the plain blue line is the regression line of the linear model of the relationship between observed and predicted stand daily ETR, with confidence interval represented with the grey dashed lines; the dashed red line is the 1:1 line. See Table S11 associated statistics. Color code for the seasons as follows:</figcaption></figure>

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

Winter;

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

Spring;

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

Summer;

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

Autumn
<figure><img src="image_178.png" /></figure>

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

<figure><img src="image_180.png" /><figcaption>**Figure W7 | Predicted versus observed soil water quantity (SWC).** For each simulation site, the plain blue line is the regression line of the linear model of the relationship between observed and predicted SWC, with confidence interval represented with the grey dashed lines; the dashed red line is the 1:1 line. See Table S12 for associated statistics. Color code for the seasons as follows:</figcaption></figure>

●, Winter; ●, Spring; ●, Summer; ●, Autumn
<figure><img src="image_181.png" /></figure>

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

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

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

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

**Figure W8** | Predicted versus observed evolution of aggregate daily species transpiration. For each simulation site, the blue line depicts the aggregated daily transpiration (mm) generated by PHOREAU from all the trees of the given species. The red line depicts the observed daily transpiration value for this species, upscaled from sapflow measurements made for individual trees using stand LAI and species leaf area to sapwood area ratios.
<figure><img src="image_186.png" /></figure>

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

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

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

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

**Figure W9 | Predicted versus observed species aggregate daily transpirations.** For each simulation site, the plain blue line is the regression line of the linear model of the relationship between observed and predicted species aggregate daily transpiration (mm), with confidence interval represented with the grey dashed lines; the dashed red line is the 1:1 line. See Table S8 for associated statistics. Color code for the seasons as follows: ●Winter; ●Spring; ●,Summer; ●,Autumn
<figure><img src="image_191.png" /></figure>

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

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

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

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

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

2060

**Figure W10 | Predicted versus observed daily stem water potential.** For each dominant species of the four simulation sites, each point represents a day with water potential observations (mPa), plotted against its corresponding predicted value by the PHOREAU model. For Puéchabon, Font Blanche and Barbeau minimum daily water potential is plotted, while the predawn potentials are shown for the Hesse site. The plain blue line is the regression line of the linear model of the relationship between observed and predicted water potential, with confidence interval represented with the grey dashed lines; the dashed red line is the 1:1 line. See Table S10 for associated statistics. Color code for the seasons as follows:
* Winter: ●,Spring: ●,Summer: ●,Autumn

2065

2070

2075

2080
<figure><img src="image_197.png" /></figure>

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

Figure W11 | Regressions of PHOREAU stand basal area increment (BAI) prediction errors against site characteristics. For each of the 340 simulation inventories, the mean annual BAI prediction error, color coded by site dominant species (see legend in bottom left) is plotted against : (a) Mean annual degree day sum (°C); (b) mean annual temperature (°C); (c) mean annual precipitation sum (mm/year); (d) soil maximum available water capacity (cm); (e) length of simulation; (f) initial inventory basal area (m²/ha). (g) initial stand density (trees/ha). The plain blue line is the regression line of the linear model of the relationship between prediction error and stand characteristic, with slope value and significance in top left, along with Pearson
<figure><img src="image_199.png" /></figure>

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

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

**Figure W12 | Predicted (by PHOREAU) against observed litter leaf area index (LAI) for available RENECOFOR inventories.** The y-axis shows the LAI predicted by the model from the stand inventory at the start of the simulation, while the x-axis represents the LAI value inferred from litter collection for the machining coordinate and closest available year. Stand points are color coded by dominant species (see legend in bottom left). The size of points shows inventory basal area. The dashed red-line is the 1:1 line; the black full line represent the regression line of the linear model between observed and predicted LAI, with confidence interval represented by the grey shaded area. Associated statistics in Table S7.
<figure><img src="image_202.png" /></figure>

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

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

**Figure W13 | Predicted (by ForCEEPS) against observed mean annual tree basal increments (BAI) for** all simulated trees over the 340 RENECOFOR and ICP II validation inventories. Tree points are color coded by species (see legend above). The dashed red-line is the 1:1 line; other full lines represent the regression lines of the linear model between observed and predicted tree productivity, with confidence intervals represented by the grey shaded area (in black the overall regression; coloured lines for species-specific regressions). Species-specific regressions are only shown for stand dominant species (in bold in legend). Associated statistics for the global simulation in top left, while species-specific statistics can be found in Table S1.
<figure><img src="image_205.png" /></figure>

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

<figure><img src="image_207.png" /><figcaption>**Figure W14 | Predicted (by ForCEEPS) against observed mean annual stand basal increments (BAI)**
for all 340 RENECOFOR and ICP II validation inventories. Stand points are color coded by dominant species (see legend above). The dashed red-line is the 1:1 line; other full lines represent the regression lines of the linear model between observed and predicted stand productivity, with confidence intervals represented by the grey shaded area (in black the overall regression; coloured lines for species-specific regressions). Associated statistics for the global simulation in top left, while species-specific statistics can be found in Table S2.</figcaption></figure>

2095