# HETEROFOR 1.0: a spatially explicit model for exploring the response of structurally complex forests to uncertain future conditions – Part 1: Carbon fluxes and tree dimensional growth

Mathieu Jonard, Frédéric André, François de Coligny, Louis de Wergifosse, Nicolas Beudez, Hendrik Davi, Gauthier Ligot, Quentin Ponette, Caroline Vincke

► To cite this version:

Mathieu Jonard, Frédéric André, François de Coligny, Louis de Wergifosse, Nicolas Beudez, et al. HETEROFOR 1.0: a spatially explicit model for exploring the response of structurally complex forests to uncertain future conditions – Part 1: Carbon fluxes and tree dimensional growth. Geoscientific Model Development, 2020, 13 (3), pp.905-935. 10.5194/gmd-13-905-2020. hal-05077035

HAL Id: hal-05077035

https://hal.science/hal-05077035v1

Submitted on 21 May 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>

# HETEROFOR 1.0: a spatially explicit model for exploring the response of structurally complex forests to uncertain future conditions – Part 1: Carbon fluxes and tree dimensional growth

Mathieu Jonard<fnref n="1" />, Frédéric André<fnref n="1" />, François de Coligny<fnref n="2" />, Louis de Wergifosse<fnref n="1" />, Nicolas Beudez<fnref n="2" />, Hendrik Davi<fnref n="3" />, Gauthier Ligot<fnref n="4" />, Quentin Ponette<fnref n="1" />, and Caroline Vincke<fnref n="1" />

**Correspondence:** Mathieu Jonard (mathieu.jonard@uclouvain.be)

Received: 12 April 2019 – Discussion started: 12 June 2019

Revised: 20 January 2020 – Accepted: 27 January 2020 – Published: 5 March 2020

**Abstract.** Given the multiple abiotic and biotic stressors resulting from global changes, management systems and practices must be adapted in order to maintain and reinforce the resilience of forests. Among others, the transformation of monocultures into uneven-aged and mixed stands is an avenue to improve forest resilience. To explore the forest response to these new silvicultural practices under a changing environment, one needs models combining a process-based approach with a detailed spatial representation, which is quite rare.

We therefore decided to develop our own model (HETEROFOR for HETEROgeneous FORest) according to a spatially explicit approach, describing individual tree growth based on resource sharing (light, water and nutrients). HETEROFOR was progressively elaborated within Capsis (Computer-Aided Projection for Strategies in Silviculture), a collaborative modelling platform devoted to tree growth and stand dynamics.

This paper describes the carbon-related processes of HETEROFOR (photosynthesis, respiration, carbon allocation and tree dimensional growth) and evaluates the model performances for three broadleaved stands with different species compositions (Wallonia, Belgium). This first evaluation showed that HETEROFOR predicts well individual radial growth (Pearson's correlation of 0.83 and 0.63 for the European beech and sessile oak, respectively) and is able to repro-

duce size-growth relationships. We also noticed that the net to gross primary production (npp to gpp) ratio option for describing maintenance respiration provides better results than the temperature-dependent routine, while the process-based (Farquhar model) and empirical (radiation use efficiency) approaches perform similarly for photosynthesis. To illustrate how the model can be used to predict climate change impacts on forest ecosystems, we simulated the growth dynamics of the mixed stand driven by three IPCC climate scenarios. According to these simulations, the tree growth trends will be governed by the CO₂ fertilization effect, with the increase in vegetation period length and the increase in water stress also playing a role but offsetting each other.

## 1 Introduction

Forest structure and composition result from soil and climate conditions, management, and natural disturbances. All of these drivers of forest ecosystem functioning are rapidly evolving due to global changes (Aber et al., 2001; Lindner et al., 2010; Campioli et al., 2012). While environmental and societal changes are taking place and will continue to happen in the future, their magnitude and the way they will occur locally remain largely uncertain (Lindner et al., 2014). Designing silvicultural systems and selecting tree species

<footnote n="1">¹Earth and Life Institute, Université catholique de Louvain, Louvain-la-Neuve, 1348, Belgium</footnote>

<footnote n="2">²AMAP, Univ Montpellier, CIRAD, CNRS, INRAE, IRD, 34000 Montpellier, France</footnote>

<footnote n="3">³Ecologie des Forêts Méditerranéennes (URFM), Institut National de la Recherche pour l'Agriculture, l'Alimentation et l'Environnement (INRAE), Avignon, 84914, France</footnote>

<footnote n="4">⁴Gembloux Agro-Bio Tech, Université de Liège, Gembloux, 5030, Belgium</footnote>
adapted to future conditions seem therefore to be risky bets (Ennos et al., 2019). Messier et al. (2015) proposed another vision of the forests in which they are considered as complex adaptive systems whose future dynamics are inherently uncertain. To maintain the ability of forests to provide a large range of goods and services in whatever future conditions, their resilience and adaptability must be improved by favouring an uneven-aged structure and tree species mixture (Thompson et al., 2009; Oliver et al., 2015). As the combinations of site conditions, climate projections, stand structures and tree species compositions are nearly infinite, all of the management options that could potentially enhance the resilience and adaptive capacity of forests cannot be tested in situ (Cantarello et al., 2017). Furthermore, such silvicultural trials provide results only for the long run, given the life span of trees, and cannot anticipate future conditions. Scenario analyses based on model simulations are therefore useful for selecting the most promising management strategies and evaluating their long-term sustainability. To explore the forest response to new silvicultural practices and yet unexperienced climate conditions in a realistic way, one needs new process-based models that are able to deal with mixed and structurally complex stands and to incorporate uncertainties in future conditions (Berger et al., 2008; Bravo et al., 2019).

In connection with the traditional forestry, viewing forests as stable systems that can be controlled, many empirical models were developed to predict tree growth in monocultures, considering that past conditions will remain unchanged in the future. On the other hand, scientists developed process-based ecophysiological models to better understand the short- and long-term forest ecosystem responses to multiple and interacting environmental changes (Dufrêne et al., 2005). This can indeed not be done through direct experimentation because the multisite and multifactorial experiments required for doing so would be too complex and too expensive (Aber et al., 2001; Boisvenue and Running, 2006). Most experiments of environment manipulation focus on single or few factors during a limited period of time, which precludes to properly take into account interactions, feedbacks and acclimation. To simplify the mathematical formalization of ecophysiological processes (e.g. radiation interception) and limit the calculation time, these process-based models were first designed for pure even-aged stands without considering the spatial heterogeneity of stand structure.

With the increasing interest in uneven-aged stands and tree species mixtures, cohort and tree-level models were also developed. Pretzsch et al. (2015) reviewed 54 forest growth models to show how they represent species mixing. Among those models, 36 were process-based with nine at the stand, 11 at the cohort and 16 at the tree level. While cohort models allow a description of the vertical structure of the stand, tree-level models are generally necessary to consider the spatial heterogeneity in the horizontal dimensions. To represent stand structure in three dimensions, the model must not only operate at the individual level but also consider the tree po-

sition. In the review of Pretzsch et al. (2015), 11 process-based models were individual-based and spatially explicit, but only three of them accounted simultaneously for radiation transfer, water cycling and phenology (i.e. BALANCE, EMILION and MAESPA). Since it describes canopy and water balance processes using a state-of-the-art approach (based on a fine crown discretization), MAESPA is a very useful tool for analysing outcomes of ecophysiological experiments (Duursma and Medlyn, 2012). MAESPA is however not suitable for multiyear simulations since it contains no routine for carbon allocation, respiration and tree dimensional growth. EMILION is also restricted to 1-year simulation (no organ emergence) and is specific to pine species with a quite detailed structural approach (Bosc et al., 2000). In contrast, tree dimensional growth is well described in BALANCE, which possesses a fine representation of tree structure (Grote and Pretzsch, 2002). In BALANCE, radiation interception by trees and water cycling are based on simpler ecophysiological concepts compared to MAESPA, and photosynthesis is calculated with a 10 d time step using the routine of Haxeltine and Prentice (1996). As the Forest v5.1 model (Schwalm and Ek, 2004), BALANCE has the advantage of merging two traditions, conventional growth and yield models together with process-based approaches, providing outputs familiar to foresters (classical tree and stand measurements obtained from forest inventory) as well as carbon fluxes and stocks. Among the three models, BALANCE is the only one that considers mineral nutrition through the impact of nitrogen (N) availability on tree growth. Some soil chemistry processes (e.g. ion exchange or mineral weathering) are however not described although they are essential to estimate bioavailability of the major nutrients other than N (P, K, Mg and Ca). Not considered in the review of Pretzsch et al. (2015), iLand is another individual-based model that describes the ecophysiological processes with an intermediate level of detail using simplified ecophysiological concepts (such as the radiation use efficiency approach) in order to also simulate forest dynamics at the landscape scale. Later, Simioni et al. (2016) developed the NOTG 3-D model to study water and carbon fluxes in Mediterranean forests using an individual-based approach to account for the spatial structure of the stand. This model is more suited for short-term simulations (a few years) than long-term (a rotation) simulations since tree dimensions are updated based on fixed empirical relationships between diameter at breast height (dbh) and tree height or crown radius.

As the models accounting for both the functional and spatial complexity are rare, we developed a new model (HETEROFOR) using a spatially explicit approach to describe individual tree growth based on resource use (light, water and nutrients) in heterogeneous forests. While the BALANCE and iLand models existed and responded roughly to our expectations, we decided to build a new model for several reasons. First, we thought that another model of this particular type would not be redundant if it was based on other con-
cepts. Instead of calculating an index of light availability, we
chose to estimate radiation interception for all trees using a
ray tracing approach. For calculating photosynthesis and tree
transpiration, we selected the Farquhar model with a shorter
time step than in BALANCE in order to account for hourly
variations in climate and soil water conditions. While we
used a slightly more complex approach for the water balance
module (Darcy approach instead of bucket model for soil wa-
ter dynamics and rainfall partitioning when passing through
the canopy), our model rests on a simpler representation of
tree structure than BALANCE. Second, we aimed at incor-
porating a detailed tree nutrition and nutrient cycling module
since we realized the necessity to integrate nutritional con-
straints in forest growth modelling, especially for predicting
the response to climate change (Fernandez-Martinez et al.,
2014; Jonard et al., 2015). Finally, we wanted to develop the
model within the frame of a collaborative modelling platform
dedicated to tree growth and stand dynamics. Among the var-
ious platforms, Capsis was the only one allowing multimodel
integration and providing a user-friendly interface (Dufour-
Kowalski et al., 2012). HETEROFOR was therefore progres-
sively elaborated through the integration of various modules
(light interception, phenology, water cycling, photosynthesis
and respiration, carbon allocation, mineral nutrition, and nu-
trient cycling) within Capsis. The advantage of such a plat-
form is to use common development environment, model ex-
ecution system, user interface and visualization tools and to
share data structures, objects, methods and libraries.

To simulate the response of forests to management and
changing environmental conditions, structuring and integrat-
ing the existing knowledge into the process-based models is
essential but not sufficient. These models must also be docu-
mented and evaluated in order to know exactly their strengths
and limits when analysing their outputs. The objectives of
this paper are to (i) describe the carbon-related processes
of HETEROFOR (photosynthesis, respiration, carbon allo-
cation and tree dimensional growth), (ii) evaluate the model
ability in reconstructing tree growth in three broadleaved
stands of different species composition and compare various
options for describing photosynthesis, respiration and crown
extension, and (iii) illustrate its potentialities by simulating
tree growth dynamics in an oak and beech stand under vari-
ous IPPC (Intergovernmental Panel on Climate Change) cli-
mate scenarios.

## 2 Materials and methods

### 2.1 Overall operation of the HETEROFOR model

HETEROFOR is a model integrated in the Capsis
(Computer-Aided Projection for Strategies in Silviculture)
platform, which is dedicated to forest growth and dynam-
ics modelling (Dufour-Kowalski et al., 2012). HETEROFOR
uses the Capsis execution system and its methods to run sim-

ulations and display the results. When running simulations
with HETEROFOR, Capsis creates a new project in which
the variables that describe the forest state are stored at a
yearly time step, starting from the initial forest characteris-
tics (initial step). Some variables (foliage state, water fluxes,
npp and gpp) are stored at an hourly or daily time step in
java objects created annually. This information is accessible
to the user through exports (see user manual in the Supple-
ment). Though some data structures and methods are shared
with other models integrated in Capsis, the initialization and
evolution procedures are specific to HETEROFOR.

For the initialization, HETEROFOR loads a series of files
containing tree species parameters, input data on trees ( loca-
tion, dimensions and chemistry), soil (chemical and physical
properties) and open-field hourly meteorological data. These
data are used to create trees and soil horizons at the initial
step. The tree is divided in three structural compartments
 branch, stem and root) and three functional ones (leaf, fine
 root and fruit). Then, HETEROFOR predicts tree growth at a
yearly time step based on underlying processes modelled at
finer time steps and at different spatial levels.

After the initialization step and at the end of each succes-
sive yearly time step, the phenological periods for each de-
ciduous species (leaf development, leaf colouring and shed-
ding) are defined for the next step from meteorological data.
When no hourly meteorological measurements are available,
the vegetation period is defined by the user who provides the
budburst and the leaf shedding dates. Knowing the key phe-
nological dates and the rates of leaf expansion, colouring and
falling, the foliage state of the deciduous species is predicted
with a daily time step during the year (de Wergifosse et al.,
2019). It is characterized by the proportions of leaf biomass
and of green leaves relative to complete leaf development,
which are key variables in the simulation of energy, water
and carbon fluxes within the forest ecosystem. The propor-
tion of green leaves impacts photosynthesis, leaf respiration
and tree transpiration, as these processes are not active any-
more on discoloured leaves, which however still intercept so-
lar radiation and rainfall. Based on a ray tracing approach,
the SamsaraLight library of Capsis (Courbaud et al., 2003)
calculates the proportions of solar radiation absorbed by the
trunk and the crown of each individual tree and the radiation
transmitted to the ground on average over the whole vegeta-
tion period (simplified radiation budget) or hourly for several
key dates (detailed radiation budget). Predicting how solar
energy is distributed within the forest ecosystem is necessary
to estimate foliage, bark and soil evaporation, tree transpira-
tion, and leaf photosynthesis.

Every hour, HETEROFOR performs a water balance and
updates the water content of each horizon. Rainfall is parti-
tioned into throughfall, stemflow and interception (André et
al., 2008a, b, 2011). Part of the rainfall directly reaches the
ground (throughfall), while the rest is intercepted by foliage
and bark. They both have a certain water storage capacity
which is regenerated by evaporation. When the foliage is sat-
urated, the overflow joins the throughfall flux whose propor-
tion increases. As the bark saturates, water flows along the
trunk to form stemflow. Throughfall and stemflow supply the
first soil horizon (forest floor) with water, while soil evapo-
ration and root uptake deplete it. The water evaporation from
the soil (as well as from the foliage and the bark) is calcu-
lated with the Penman–Monteith equation based on the solar
radiation absorbed by each component. Using the same equa-
tion, individual tree transpiration is estimated by determining
the stomatal conductance from tree characteristics, soil wa-
ter potential and meteorological conditions. The distribution
of root water uptake among the soil horizons is done accord-
ing to the soil water potential and the vertical distribution of
fine roots. Water exchanges between soil horizons are con-
sidered as water inputs (capillary rise) or outputs (drainage).
These soil water transfers are calculated based on the soil
water potential gradient according to the Darcy law and us-
ing pedotransfer functions to determined soil hydraulic prop-
erties. By default, HETEROFOR calculates the water fluxes
at the stand scale by aggregating individual fluxes (i.e. tree
transpiration) or tree properties (e.g. foliage and bark capac-
ity or stemflow proportion). With this option, all trees take
up water in the same soil horizons, assuming that soil wa-
ter is redistributed homogeneously between two hourly time
steps. However, the user can choose an alternative option to
calculate all of the water fluxes at the individual level. In this
case, the model distributes the total soil volume in individ-
ual soil volumes (called pedons) and performs a water bal-
ance for each one. Contrary to the default option, assuming a
homogeneous horizontal water redistribution, the alternative
option supposes no water redistribution among pedons (de
Wergifosse et al., 2019).

The user can choose to calculate the gross primary pro-
duction of each tree (gpp) either based on a radiation use
efficiency approach distinguishing sunlit and shaded leaves
(yearly time step) or using the Farquhar et al. (1980) model
(hourly time step). The latter is analytically coupled to the
stomatal conductance model proposed by Ball et al. (1987).
The photosynthesis is computed using the library CAS-
TANEA that is also present in Capsis (Dufrêne et al., 2005).
This calculation requires the proportions of sunlit and shaded
leaves, the direct and diffuse photosynthetically active radia-
tion (PAR) absorbed per unit leaf area, and the mean soil wa-
ter potential. At the end of the vegetation period, gpp is con-
verted to net primary production (npp) after the subtraction
of growth and maintenance respiration. Maintenance respi-
ration is either considered as a proportion of gpp (depending
on the crown to stem diameter ratio) or calculated hourly for
each tree compartment by considering the living biomass,
the nitrogen concentration and a Q₁₀ function for the tem-
perature dependency following Ryan (1991) as in Dufrêne et
al. (2005). Carbon allocation is done once a year at the end
of the vegetation period, which allows an update of the tree
dimensions for the next yearly time step, during which tree
size does not change. Allocation of carbon to foliage and fine

roots is prioritized and carried out by ensuring a functional
balance between carbon fixation and nutrient uptake through
a fine root to leaf biomass ratio that depends on the tree nu-
tritional status (Helmisaari et al., 2007). Allometric relation-
ships are then used to describe carbon allocation to structural
components (trunk, branches and structural roots) and to de-
rive tree dimensional growth (diameter at breast height, total
height, height to crown base, height of maximum crown ex-
tension and crown radii in four directions) while considering
competition with neighbouring trees (Fig. 1).

Knowing the chemical composition of the tree compart-
ments for a given tree nutrient status, HETEROFOR com-
putes the individual tree nutrient requirements based on the
estimated annual growth rate and deduces the tree nutrient
demand after a subtraction of the amount of retranslocated
nutrients. In parallel, the potential nutrient uptake (soil nutri-
ent supply) is obtained by calculating the maximum rate of
ion transport towards the roots (by diffusion and mass flow).
The actual uptake is then determined by adjusting the tree
nutrient status and growth rate so that tree nutrient demand
matches soil nutrient supply. The nutrient limitation of tree
growth is achieved through the regulation of photosynthesis,
maintenance respiration and through the effect of the tree nu-
trient status on fine root allocation.

The soil chemistry is characterized at the tree or stand
scale for the various soil horizons defined by the user. In each
soil horizon, the chemical composition of the soil solution is
in equilibrium with the exchange complex and the secondary
minerals. Soil receives nutrients coming from atmospheric
deposition, organic matter mineralization and primary min-
eral weathering and is depleted by root uptake and immobi-
lization in microorganisms. The chemical equilibrium within
the soil solution, with the exchange complex or the minerals,
is updated yearly with the PHREEQC geochemical model
(Charlton and Parkhurst, 2011) coupled to HETEROFOR
through a dynamic link library.

In this paper, we present a detailed description of the pro-
cesses regulating the carbon fluxes (Fig. 1), while the phenol-
ogy and water balance modules are presented in a companion
paper (de Wergifosse et al., 2019) and the nutrient cycling
and tree nutrition module will be described later in a third
paper.

## 2.2 Detailed model description

### 2.2.1 Initialization

To initialize HETEROFOR, the relative position (x, y, z) and
the main dimensions of each tree must be provided, includ-
ing the following: girth at breast height (gbh; in centime-
tres), height (h; in metres), height of maximum crown ex-
tension (hlce; in metres), height to crown base (hcb; in me-
tres) and crown radii in the four cardinal directions (cr; in
metres). During the initialization phase, the biomass of each
tree compartment is calculated according to the equations
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Conceptual diagram of the HETEROFOR model. The incident PAR radiation is absorbed by individual trees using a ray-tracing model (SamsaraLight library). Then, the absorbed PAR (aPAR) is converted into gross primary production (gpp) based on the PAR use efficiency concept (first option) or with a biochemical model of photosynthesis (second option). The photosynthesis calculation depends on the soil water potential, which is updated hourly thanks to the water balance module described in detail in de Wergifosse et al. (2019). The net primary production (npp) is obtained using a npp to gpp ratio or by subtracting the growth and maintenance respiration (the latter being temperature dependent). npp is first allocated to foliage using an allometric equation function of tree diameter (dbh) and crown radius (cr). All of these processes (radiation interception, photosynthesis, respiration and evapotranspiration) depend on the foliage development stage, which is determined based on the phenology module. The carbon allocated to fine roots is determined based on a fine root-to-foliage ratio, dependent on the tree nutritional status. Fruit production is calculated with an allometric equation based on dbh and on light availability. The remaining carbon is allocated to structural compartments (roots, trunk and branches) using a fixed proportion for the below-ground part. dbh and height growth (Δdbh and Δh) are deduced from the change in aboveground biomass by deriving and rearranging an allometric equation. Finally, crown extension is predicted with a distance-dependent or distance-independent approach.</figcaption></figure>

used for carbon allocation (see Sect. 2.2.4. If available, site-specific allometric equations can also be used to calculate initial biomasses of tree compartments. When data on fruit litterfall are available, a file providing the amount of fruit litterfall per year and per tree species can be loaded and used to adapt the allometric equations predicting fruit production at the individual level. When the water balance module is activated, two additional files must be loaded, including a file describing soil horizon properties and another one for the hourly meteorology. Finally, the user must provide the nutrient concentrations of the current leaves (N, P, K, Ca and Mg) for each tree species. These foliar concentrations are then used to estimate the tree nutrient status for each major nutrient. When the tree nutrition and nutrient cycling module is not activated, these concentrations are kept constant throughout the simulation.

### 2.2.2 Gross primary production

The annual gross primary production of each tree (gpp; in kilograms of C per year) is calculated either based on a PAR use efficiency (PUE) approach (Monteith, 1977) or using the photosynthesis method of the CASTANEA model (Dufrêne et al., 2005). For the first option, the only input needed by the model is the mean monthly global radiation. The second

option requires hourly meteorological data and the activation of the water balance calculation. In any case, a series of intermediate variables are needed to calculate gpp.

For the PUE approach, the model uses the solar radiation absorbed by each tree during the vegetation period (aRAD; in megajoules per year); aRAD is then converted in PAR (aPAR; in moles of photons per year) by supposing that 46% of the solar radiation (RAD) is PAR and 1 MJ is equivalent to 4.55 moles of photons. The diffuse and direct components of aPAR are also considered (aPAR<sub>diff</sub> and aPAR<sub>dir</sub>; in moles of photons per year). While all of the leaves receive diffuse PAR, only sunlit leaves absorb direct PAR. To estimate the sunlit-leaf proportion (Prop<sub>sl</sub>) at the tree level, HETEROFOR uses an adaptation of the classical stand-scale approach based on the Beer-Lambert law (Teh, 2006):

$$ \text{Prop}_{\text{sl}} = \frac{1 - \exp(-k \cdot \text{LAI})}{k}, \quad (1) $$

with k, the extinction coefficient, and LAI, the leaf area index (in square metres per square metre).

At the individual scale, the leaf area index is calculated by dividing the tree leaf area (a<sub>leaf</sub>; in square metres) by the crown projection area (cpa; in square metres). The value obtained is then multiplied by the light competition index (LCI; in megajoules per megajoule) to account for the shading ef-
fect of the neighbouring trees as follows:

$$ \text{Prop}_{\text{sl}} = \frac{1 - \exp\left(-k \cdot \frac{a_{\text{leaf}}}{\text{cpa}}\right)}{k} \cdot \text{LCI}, \qquad (2) $$

where LCI is the ratio between the absorbed radiation calculated with and without neighbouring trees in SamsaraLight. LCI ranges from 0 (no light reaching the tree) to 1 (no light competition).

To adapt the PAR use efficiency (PUE) concept at the tree level, we considered a distinct PUE for sunlit (sl) and shaded (sh) leaves and calculated an average PUE weighted as follows:

$$ \text{pue} = \frac{a\text{PAR}_{\text{diff}} \cdot (\text{Prop}_{\text{sl}} \cdot \text{PUE}_{\text{sl}} + \text{Prop}_{\text{sh}} \cdot \text{PUE}_{\text{sh}}) + a\text{PAR}_{\text{dir}} \cdot \text{PUE}_{\text{sl}}}{a\text{PAR}} \qquad (3) $$

This pue value is then used to calculate gpp based on aPAR and a reducer accounting for water stress (red<sub>water</sub>) as follows:

$$ \text{gpp} = a\text{PAR} \cdot \text{pue} \cdot \text{red}_{\text{water}} \qquad (4) $$

The default value of red<sub>water</sub> is 1, but, when the water balance module is activated, it is set to the ratio between the actual and the potential (i.e. considering no soil water limitation) tree transpiration (t<sub>actual</sub> and t<sub>pot</sub>; in litres per year). This ratio estimates the fraction of the vegetation period during which stomata are partially or totally closed due to a limitation in soil water availability. Since this ratio is always lower or equal to 1, a correction factor is applied to avoid introducing a bias.

$$ \text{red}_{\text{water}} = \frac{t_{\text{actual}}}{t_{\text{pot}}} \cdot \text{corr.} \qquad (5) $$

gpp can also be estimated using the photosynthesis method of CASTANEA (Dufrêne et al., 2005). This method consists in the biochemical model of Farquhar et al. (1980) analytically coupled with the approach of Ball et al. (1987) that linearly relates stomatal conductance to the product of the carbon assimilation rate by the relative humidity. The slope of this relationship varies between 0 and 1, with the soil water availability characterized in HETEROFOR based on a decreasing exponential function of the mean soil water potential (see eq. 55 in de Wergifosse et al., 2019). The formulation of Ball et al. (1987) was slightly adapted to the tree level by accounting for the influence of tree height. Indeed, leaf water potential increases with leaf height and induces a decrease in stomatal conductance (Ryan and Yoder, 1997; Schäfer et al., 2000). In eq. (55) in de Wergifosse et al. (2019), stomatal conductance is inversely proportional to the height of maximum crown extension.

The photosynthesis routine requires, at an hourly time step, the direct and diffuse PAR absorbed per unit leaf area.

The direct PAR is intercepted only by sunlit leaves and is obtained by multiplying the hourly incident PAR (micromoles of photons per square metre per second) by the proportion of direct PAR absorbed by sunlit leaves. For a tree, this proportion is fixed by default for the whole vegetation period and calculated as the ratio between the direct PAR absorbed per unit sunlit-leaf area during the vegetation period (in moles of photons per square metre per year) and the incident PAR accumulated over the same period (in moles of photons per square metre per year). A similar procedure is used for the diffuse absorbed PAR, except that it is related to the total leaf area. When using the detailed version of SamsaraLight, the proportions of direct/diffuse PAR absorbed per unit leaf area change every hour during the day and depending on the physiological stage. The photosynthesis routine of CASTANEA also requires the foliar nitrogen concentration to estimate the maximal carboxylation rate (Dufrêne et al., 2005).

### 2.2.3 Growth and maintenance respiration

The value of gpp is converted to annual net primary production (npp; in kilograms of C per year) using either a ratio depending on the crown to stem diameter ratio (Eq. 6) or, after subtraction of growth (gr) and maintenance respiration (mr) (Eq. 7), according to the theory of respiration developed by Penning de Vries (1975).

$$ \text{npp} = \text{gpp} \cdot r_{\text{npp\_gpp}}(\text{DdIndex}) \qquad (6) $$

$$ \text{npp} = \text{gpp} - \text{mr} - \text{gr} \qquad (7) $$

Mäkelä and Valentine (2001) showed that the npp to gpp ratio changes with some tree characteristics (tree height and age). Based on simulated gpp and npp reconstructed by using the model in reverse mode (see Sect. 2.2.7), we tested the impact of several variables characterizing tree dimensions and shape (height, dbh, crown radius, crown volume, crown to stem diameter ratio, and aboveground volume or biomass) on the npp to gpp ratio. The best relationship was obtained with the crown to stem diameter ratio (Dd; in metres per metre), which had a negative effect on the npp to gpp ratio. This indicates that the proportion of gpp lost by respiration increases for trees with a large crown. Unfortunately, the crown to stem diameter ratio not only varies with the tree shape reflecting past competition conditions but also changes during the course of the tree development for some tree species. Therefore, we standardized it to remove the size effect in order to obtain an index (DdIndex) that only characterizes the tree shape. This index is particularly useful in accounting for the large differences in the oak crown extensions according to the silvicultural system (large crowns in former coppices with standards vs. narrow crowns in dense high forests).

$$ r_{\text{npp\_gpp}} = \alpha + \beta \cdot \text{DdIndex}, \qquad (8) $$

where α and β are parameters and DdIndex is defined as

$$ \text{DdIndex} = \frac{\text{Dd}}{\text{Dd}_{\text{pred}}}, \qquad (9) $$
with Dd, the crown to stem diameter ratio determined from the tree mean crown radius (cr<sub>mean</sub>; in metres) and diameter at breast height (dbh; in metres), and Dd<sub>pred</sub>, the crown to stem diameter ratio predicted based on the girth at breast height (gbh; in centimetres):

$$Dd_{pred} = \alpha + \beta \cdot gbh + \gamma \cdot \frac{1}{gbh} + \delta \cdot \frac{1}{gbh^2}. \quad (10)$$

In Eq. (7), maintenance respiration is calculated for each tree by summing the maintenance respiration of each compartment, which is estimated from the nitrogen content of its living biomass, and considering a Q<sub>10</sub> function for the temperature dependency. During daytime, the inhibition of foliar respiration by light is taken into account by considering that this inhibition reduces respiration by 62% (Villar et al., 1995).

$$mr = \sum_{comp.} \left( b_{comp.} \cdot f_{living} \cdot [N] \cdot R_{T_{ref}} \cdot Q_{10\_comp.}^{\frac{T-T_{ref}}{10}} \right), \quad (11)$$

with b<sub>comp.</sub>, the tree compartment biomass (kilograms of organic matter), f<sub>living</sub>, the fraction of living biomass, [N], the nitrogen concentration (in grams per kilogram), R<sub>T<sub>ref</sub></sub>, the maintenance respiration per gram of N at the reference temperature (15 °C), and T, the air temperature for aboveground tree compartments or the soil horizon temperature for roots (see Appendix A). Root maintenance respiration is estimated for each soil horizon separately.

The fraction of living biomass is fixed to 1 for leaves and fine roots or equals the proportion of sapwood for the structural tree compartments. The sapwood proportion is derived from the sapwood area (a<sub>sapwood</sub>; in square centimetres), which is determined based on an empirical function of the tree compartment diameter (Ø<sub>comp.</sub>; in centimetres) as follows:

$$a_{sapwood} = a + b \cdot \emptyset_{comp.} + c \cdot \emptyset_{comp.}^2. \quad (12)$$

Growth respiration (gr) is the sum of the tree compartment growth respirations, which are proportional to their biomass increments (see Sect. 2.2.4).

$$gr = \sum_{comp.} (R_{gr} \cdot \Delta b_{comp.}), \quad (13)$$

where R<sub>gr</sub> is the growth respiration per unit biomass increment (in kilograms of C per kilogram of C).

### 2.2.4 Carbon allocation and dimensional growth

For each tree, the npp and the carbon retranslocated from leaves and roots (rt<sub>leaf</sub> and rt<sub>fine root</sub>; in kilograms of C per year) are distributed among the various tree compartments at the end of the year. rt<sub>leaf</sub> and rt<sub>fine root</sub> are determined as follows:

$$rt_{leaf} \text{ or fine root} = b_{leaf} \text{ or fine root} \cdot \delta_{leaf} \text{ or fine root} \cdot rtr_{leaf} \text{ or fine root}, \quad (14)$$

where b<sub>leaf</sub> and b<sub>fine root</sub> are the tree leaf and fine root biomasses (in kilograms of C), δ<sub>leaf</sub> and δ<sub>fine root</sub> are the leaf and fine root turnover rates (in kilograms of C per kilogram of C per year), and rtr<sub>leaf</sub> and rtr<sub>fine root</sub> are the leaf and fine root retranslocation rates (in kilograms of C per kilogram of C).

b<sub>leaf</sub> is estimated with an allometric equation based on the stem diameter at breast height (dbh; in centimetres) and on the crown to stem diameter ratio (Dd),

$$b_{leaf} = \alpha \cdot dbh^{\beta} \cdot Dd^{\gamma}. \quad (15)$$

b<sub>fine root</sub> is deduced from the leaf biomass using the fine root to leaf ratio (r<sub>fine root\_leaf</sub>) as follows:

$$b_{fine\ root} = b_{leaf} \cdot r_{fine\ root\_leaf}. \quad (16)$$

r<sub>fine root\_leaf</sub> takes a value between a minimum (r<sub>fine root\_leaf\_min</sub>) and maximum (r<sub>fine root\_leaf\_max</sub>) ratio, depending on the tree nutritional status, in accordance with the concept of functional balance (Mäkela, 1986). This means that a higher ratio (more carbon allocation to fine roots) is used when tree suffers from nutrient deficiency. For each nutrient, a candidate ratio is obtained based on a linear relationship depending on the nutritional status. The ratio increases when the nutritional status deteriorates, and this effect is more pronounced for nitrogen (N), with nitrogen (N) > phosphorus (P) > potassium (K) > magnesium (Mg) > calcium (Ca). Among the candidate ratios, the maximum is retained in order to account for the fact that the most limiting nutrient has the dominant effect. For each nutrient, the nutritional status is bounded between 0 and 1 and calculated based on the foliar concentrations (provided in the inventory file) and the optimum and deficiency thresholds (Mellert and Göttlein, 2012).

$$Status(Nutrient) = \frac{[\text{Foliar Nutrient}] - \text{Deficiency}}{\text{Optimum-Deficiency}}. \quad (17)$$

The leaf and fine root litter amounts (s<sub>leaf</sub> and s<sub>fine root</sub>; in kilograms of C per year) are estimated based on the turnover rate taking into account the retranslocation as follows:

$$s_{leaf \text{ or fine root}} = b_{leaf \text{ or fine root}} \cdot \delta_{leaf \text{ or fine root}} \cdot (1 - r_{leaf \text{ or fine root}}). \quad (18)$$

Allocation priority is given to leaves and fine roots. The carbon allocated to leaves corresponds to the annual leaf production (p<sub>leaf</sub>; in kilograms of C per year), which is equal to the fallen leaf biomass of the previous year plus the leaf biomass change (Δb<sub>leaf</sub>; in kilograms of C per year),

$$p_{leaf} = b_{leaf_{t-1}} \cdot \delta_{leaf} + \Delta b_{leaf}, \quad (19)$$

where Δb<sub>leaf</sub> is determined by

$$\Delta b_{leaf} = b_{leaf_{t}} - b_{leaf_{t-1}}, \quad (20)$$
where $b_{leaf_{t-1}}$ and $b_{leaf_t}$ are the tree leaf biomasses corresponding to the previous and the current years, respectively.

The fine root production is then estimated according to the same logic.

$$p_{\text{fine root}} = b_{\text{fine root}_{t-1}} \cdot \delta_{\text{fr}} + \Delta b_{\text{fr}}, \quad (21)$$

where $b_{\text{fine root}_{t-1}}$ is provided by Eq. (16).

When the carbon allocated to leaves and fine roots is higher than the npp plus the retranslocated carbon (suppressed trees with low gpp and npp for their size), the leaf and fine root productions are recalculated so that they do not exceed 90% of the available carbon.

Then, the fruit production ($p_{\text{fruit}}$; in kilograms of C per year) is estimated with an allometric equation similar to Eq. (15) and is considered directly proportional to the light competition index since fructification is known to be favoured when tree crowns are exposed to the sun (Greene et al., 2002; Davi et al., 2016). A threshold dbh ($\text{dbh}_{\text{threshold}}$; in centimetres) is fixed, below which no fruit production occurs.

$$p_{\text{fruit}} = \alpha \cdot \text{LCI} \cdot (\text{dbh} - \text{dbh}_{\text{threshold}})^{\beta} \quad (22)$$

In this equation, the parameter $\alpha$ takes a default value or is adapted based on the fruit production of the year (when the file with the amount of fruit litterfall per year and per tree species is loaded).

Part of the carbon is also used to compensate for branch and root mortality. The branch mortality ($s_{\text{branch}}$; in kilograms of C per year) is described with an equation of the same form as Eq. (15), while the structural root mortality ($s_{\text{root}}$; in kilograms of C per year) is obtained using a turnover rate similar to that of the branches.

After subtracting the leaf, fine root and fruit productions and the root and branch senescence, the remaining carbon is allocated to structural tree compartment growth:

$$\Delta b_{\text{structural}} = \text{npp} + \text{rt} - p_{\text{leaf}} - p_{\text{fine root}} - p_{\text{fruit}} - s_{\text{branch}} - s_{\text{root}}. \quad (23)$$

At this stage, the remaining carbon is partitioned between the above- and below-ground parts of the tree according to a fixed root-to-shoot ratio ($r_{\text{root\_shoot}}$) as follows:

$$\Delta b_{\text{structural\_above}} = \frac{\Delta b_{\text{structural}}}{(1 + r_{\text{root\_shoot}})} \quad (24)$$

$$\Delta b_{\text{structural\_below}} = \Delta b_{\text{structural}} - \Delta b_{\text{structural\_above}}. \quad (25)$$

The increment in aboveground structural biomass is then used to determine the combined increment in dbh and total height ($h$; in metres) based on an allometric equation used to predict aboveground woody biomass (Genet et al., 2011; Hounzandji et al., 2015),

$$b_{\text{structural\_above}} = \alpha + \beta (\text{dbh}^2 \cdot h)^{\gamma}. \quad (26)$$

Deriving this equation and rearranging terms gives

$$\Delta b_{\text{structural\_above}} = \beta \gamma (\text{dbh}^2 \cdot h)^{\gamma-1} \Delta (\text{dbh}^2 \cdot h) \quad (27)$$

$$\Delta (\text{dbh}^2 \cdot h) = \frac{\Delta b_{\text{structural\_above}}}{\beta \gamma (\text{dbh}^2 \cdot h)^{\gamma-1}}. \quad (28)$$

The development of the left term provides

$$\Delta (\text{dbh}^2 \cdot h) = (\text{dbh} + \Delta \text{dbh})^2 \cdot (h + \Delta h) - \text{dbh}^2 \cdot h, \quad (29)$$

which can be further developed (see Appendix B for details) to isolate $\Delta h$:

$$\Delta h \approx \frac{\Delta (\text{dbh}^2 \cdot h)}{\text{dbh}^2} - \frac{h \cdot \Delta \text{dbh}^2}{\text{dbh}^2} \quad (30)$$

From Eq. (30), we know that the height increment can be expressed as a function of $\frac{\Delta (\text{dbh}^2 \cdot h)}{\text{dbh}^2}$. In the following, we refer to it as the height growth potential ($\Delta h_{\text{pot}}$) since it corresponds to the height increment if all of the remaining carbon was allocated to height growth. Contrary to the other term of Eq. (30) ($\frac{h \cdot \Delta \text{dbh}^2}{\text{dbh}^2}$), which is unknown, this height growth potential can be evaluated at this step by dividing the result of Eq. (28) by $\text{dbh}^2$. However, depending on the level of competition for light and on the tree size, only part of this height growth potential will be effectively realized for the height increment. For each tree species, an empirical relationship, predicting height growth from the height growth potential, the light competition index and the tree size (dbh or height), was therefore fitted based on successive inventory data (see Appendix E).

$$\Delta h = a + b \cdot \text{dbh} + c \cdot h + d \cdot \text{LCI} + e \cdot \Delta h_{\text{pot}} + f \cdot (\Delta h_{\text{pot}})^2 + g \cdot (\Delta h_{\text{pot}})^3 \quad (31)$$

The dbh increment is then determined by rearranging Eq. (29) as follows:

$$\Delta \text{bh} = \sqrt{\frac{\Delta (\text{dbh}^2 \cdot h) + \text{dbh}^2 \cdot h}{(h + \Delta h)}} - \text{bh}. \quad (32)$$

The increments in root, stem and branch biomasses are obtained as follows:

$$\Delta b_{\text{root}} = r_{\text{root\_shoot}} \cdot \Delta b_{\text{structural\_above}} \quad (33)$$

$$\Delta b_{\text{stem}} = f \cdot \rho \cdot \left( (\text{dbh} + \Delta \text{dbh})^2 \cdot (h_{\text{del}} + \Delta h_{\text{del}}) - \text{dbh}^2 \cdot h_{\text{del}} \right) \quad (34)$$

$$\Delta b_{\text{branch}} = \Delta b_{\text{structural\_above}} - \Delta b_{\text{stem}}, \quad (35)$$

where $f$ is the form coefficient (in cubic metres per cubic metre), $\rho$ is the stem volumetric mass (in kilograms of C per cubic metre) and $h_{\text{del}}$ is the Delevoy height (in metres) corresponding to the height at which stem diameter is half the diameter at breast height (see Appendix C).
The branch and root biomasses are then distributed in three categories, defined based on the diameter, which are as follows: small branches/roots < 4 cm, medium branches/roots between 4 and 7 cm and coarse branches/roots > 7 cm. The proportions of small, medium and coarse branches/roots are determined based on equations of the same form as those presented in Hounzandji et al. (2015) for oak branches. Until we can adjust these equations on appropriate data sets, the parameters of Hounzandji et al. (2015) are also used for beech branches and for oak and beech roots. The distribution in the root categories has no impact on the functioning of the model since this information is not used elsewhere. This is just a model output that the user can ignore or consider as a whole.

### 2.2.5 Crown extension

Depending on whether the competition with the neighbouring trees is taken into account or not, the crown dynamics can be described by two different approaches. When local competition is not considered (distance-independent approach), changes in crown dimensions are derived from dbh or height increment based on the following empirical relationships:

$$ \Delta hlce = hlce\% \cdot \Delta h \qquad (36) $$

$$ \Delta hcb = hcb\% \cdot \Delta h \qquad (37) $$

$$ \Delta cr = Dd_{pred} \cdot \frac{\Delta dbh}{200}, \qquad (38) $$

where hcb% and hlce% are the proportions of the total height corresponding to the height to crown base (hcb; in metres) and to the height of largest crown extension (hlce; in metres), respectively; Δcr is the change in crown radius (in metres) whatever the direction; Dd$_{pred}$ is the crown to stem diameter ratio estimated by Eq. (10).

Alternatively, the changes in crown dimensions can be described based on the competition with the neighbouring trees (distance-dependent approach). The space around a target tree is divided into four sectors according to the four cardinal directions (north between 315 and 45°, east between 45 and 135°, south between 135 and 225°, and west between 225 and 315°). In each sector, the tree that is the closest to the target tree is retained as a competitor if its height is higher than the hcb of the target tree. Beyond a certain distance (i.e. 2 times the maximal crown radius of 10 m), no competitor is considered. For each main direction, the model calculates an hlce at equilibrium (hlce$_{eq}$; in metres) for the target tree. This hlce at equilibrium is located between a minimum (hcb; in metres) and a maximum (hlce$_{max}$; in metres). hlce$_{max}$ is obtained by determining the higher intersection between the potential crowns of the target tree and the competitor. The potential crown of a tree is the crown that this tree would have had in absence of competition and is considered as having the shape of a half-ellipsoid, centred on the tree trunk and with the semi-axis lengths equal to the tree potential crown

radius (cr$_{pot}$; in metres; see below) and to the crown length (h-hcb). hlce$_{eq}$ is positioned between the minimum and the maximum values according to the competition intensity, estimated based on the target tree and the competitor heights (h$_{target}$ and h$_{comp}$; in metres), as well as the hcb of the target tree (Appendix D).

$$ hlce_{eq} = hcb + (hlce_{max} - hcb) \cdot \max \left( 0, \min \left( 1, \frac{h_{comp} - hcb}{h_{target} - hcb} \right) \right) \qquad (39) $$

The four values of hlce$_{eq}$ are then averaged (hlce$_{eq\_mean}$).

Finally, the change in hlce is determined as follows: if hlc < hlce$_{eq\_mean}$,

$$ \Delta hlce = \min(\Delta hlce_{max}, hlce_{eq\_mean} - hlce) \qquad (40) $$

else,

$$ \Delta hlce = \max(-\Delta hlce_{max}, hlce_{eq\_mean} - hlce), \qquad (41) $$

where Δhlce$_{max}$ is the maximum change in hlc allowed by the model.

The change in hcb is obtained with the same logic. If hcb < hcb$_{eq\_mean}$,

$$ \Delta hcb = \min(\Delta hcb_{max}, hcb_{eq\_mean} - hcb), \qquad (42) $$

else,

$$ \Delta hcb = \max(-\Delta hcb_{max}, hcb_{eq\_mean} - hcb), \qquad (43) $$

where hcb$_{eq\_mean}$ is the hcb estimated from the tree height based on hcb% (Eq. 37).

The changes in the four crown radii are calculated based on crown radii at equilibrium (cr$_{eq}$; in metres), which are estimated by considering the competitive strength of the target and neighbouring trees. For a given direction, cr$_{eq}$ is calculated based on the potential (free-growth) crown radius of the target tree (cr$_{pot\_target}$; in metres) and its competitor (cr$_{pot\_comp}$; in metres), the distance between the two trees (d; in metres), and the crown overlap ratio (r$_{overlap}$; in metres per metre) as follows:

$$ cr_{eq} = \frac{cr_{pot\_target}}{cr_{pot\_target} + cr_{pot\_comp}} \cdot d \cdot r_{overlap\_target} \qquad (44) $$

The potential crown radius (cr$_{pot}$) of a tree if determined by

$$ cr_{pot} = \frac{dbh}{200} \cdot Dd_{pred} \cdot sh, \qquad (45) $$

where Dd$_{pred}$ is the crown to stem diameter ratio estimated by Eq. (10), and sh is a coefficient allowing a shift from the mean to the maximum Dd$_{pred}$.

The crown overlap ratio is estimated by considering neighbouring trees of the same species two by two and then calculating the ratio between the sum of their crown radii and the
distance between the corresponding tree stems. This overlap ratio accounts for the capacity of a tree species to penetrate in neighbouring crowns.

The change in crown radius is then determined for each direction as follows:

if $cr < cr_{eq}$,

$$ \Delta cr = \min(\Delta cr_{max}cr_{eq} - cr), \quad (46) $$

else,

$$ \Delta cr = \max(-\Delta cr_{max}cr_{eq} - cr), \quad (47) $$

where $\Delta cr_{min}$ and $\Delta cr_{max}$ are respectively the minimum and the maximum change in cr allowed by the model. They are obtained in a similar way as $cr_{pot}$, with

$$ \Delta cr_{max} = \frac{\Delta dbh}{200} \cdot Dd \cdot sh. \quad (48) $$

### 2.2.6 Tree harvesting and mortality

During the simulation, thinning can be achieved at each annual step by (i) selecting the trees from a list or a map or according to tree characteristics (tree species, age, dbh, height, etc.), (ii) defining the number of trees to be thinned per diameter class using an interactive histogram, or (iii) loading a file listing the trees that must be thinned. In addition, the thinning methods developed for GYMNOS and QUERGUS are compatible with HETEROFOR. They allow one to reach a target basal area, density or relative density index by thinning from below or from above or by creating gaps (Ligot et al., 2014).

When the npp of a tree is not sufficient to ensure a normal leaf and fine root development (for suppressed trees and/or after a severe drought), the leaf biomass is reduced and induces a defoliation, which is estimated as follows:

$$ Def = \frac{b_{leaf} - b_{leaf\_corr}}{b_{leaf}} \cdot 100, \quad (49) $$

where $b_{leaf}$ and $b_{leaf\_corr}$ are, respectively, the leaf biomass estimated with Eq. (15) and the leaf biomass corrected to match the available carbon (see Sect. 2.2.4).

Tree mortality occurs when trees reach a defoliation of 90%, considering that a tree with less than 10% of its leaves is in an advanced stage of decline and is unlikely to recover (Manion, 1981). Hence, HETEROFOR takes into account the mortality resulting from carbon starvation due to light competition and/or water stress (stomatal closure).

### 2.2.7 Growth reconstruction

HETEROFOR was adapted to allow the user to run it in reverse mode by starting from the known increments in dbh and h to reconstruct an individual npp using exactly the same parameters and equations as in the normal mode. To achieve a reconstruction, an inventory file with tree measurements

must be loaded to create the initial step. From this initial step, the reconstruction tools can be launched and require another inventory file with tree measurements achieved one or several years later. Based on these two inventories, HETEROFOR calculates the mean dbh and h increments for each tree and uses the model equations to reconstruct each step and evaluate among other individual npp's. The npp is obtained by rearranging Eq. (23), in which the carbon allocated to the structural biomass is calculated from the dbh and h increments using Eqs. (27), (25) and (24). The carbon allocated to leaf, fine root and fruit production is determined, respectively, with Eqs. (19), (21) and (22), while the amount retranslocated from leaves and roots before senescence is evaluated with Eq. (14). Finally, the terms of Eq. (23) accounting for the leaf and fine root litter were determined with Eq. (18). In addition to two stand inventories, the reconstruction tool also requires a file listing the trees which were cut or died between the two inventory dates and the last year during which they were present in the stand.

## 2.3 Input variables and parameter setting for a case study

The model was tested in three stands that contrast in structure and species composition. These stands were located close to each other (< 1 km) on the same tableland (300 m elevation) in the western part of the Belgian Ardennes at Baileux (50°01'N, 4°24'E). The average annual rainfall is slightly above 1000 mm, and the mean annual temperature is 8 °C. The forest (60 ha) consists of sessile oaks (Quercus petraea Liebl.) and European beeches (Fagus sylvatica L.) and lies on acid brown earth soil (Luvisol according to the FAO soil taxonomy) with a moder humus and an $A_h B_w C$ profile. The soil has been developed on a loamy and stony solifluction sheet in which weathering products of the bedrock (Lower Devonian: sandstone and schist) were mixed with added periglacial loess.

By the end of the 19th century, the Baileux forest was probably an oak coppice with a few standards. Taking advantage of the massive oak regeneration in the 1880s, the forest developed progressively into a high forest and was then invaded by beeches. In 2001, the area was covered by aged oak trees and heterogeneously sized beech trees. At that time, three experimental plots were installed at the Baileux site in order to study the impact of tree species mixing on ecosystem functioning (Jonard et al., 2006, 2007, 2008; André et al., 2008a, b, 2010, 2011); two plots were located in stands dominated by either sessile oak or beech, and the third one was a mixture of both species (Table 1). In each plot, all trees with a circumference higher than 15 cm were mapped (coordinates) and measured (stem circumference at a height of 1.3 m, total tree height, height of largest crown extension, height to crown base and crown diameters in two directions) at the end of the years 2001 and 2011.
<table><caption>**Table 1.** Stand characteristics for the main tree species derived from stand inventories in 2001. Standard deviation is provided in parentheses.</caption><thead><tr><th>Stand</th><th>Tree species</th><th>Tree density<br>(N ha<sup>-1</sup>)</th><th>Basal area<br>(m<sup>2</sup> ha<sup>-1</sup>)</th><th>gbh*<br>(cm)</th><th>Dominant height<br>(m)</th></tr></thead><tbody><tr><td>Oak-dominated<br>(0.90 ha)</td><td>Sessile oak</td><td>187</td><td>16.2</td><td>100.6 (26.5)</td><td>21.9</td></tr><tr><td></td><td>European beech</td><td>118</td><td>4</td><td>46.4 (35.6)</td><td>15.5</td></tr><tr><td>Beech-dominated<br>(1.44 ha)</td><td>Sessile oak</td><td>72</td><td>6.4</td><td>103.3 (18.1)</td><td>23</td></tr><tr><td></td><td>European beech</td><td>217</td><td>16.5</td><td>87.5 (41.5)</td><td>25</td></tr><tr><td>Mixed<br>(1.80 ha)</td><td>Sessile oak</td><td>118</td><td>12.9</td><td>115.5 (21.0)</td><td>24.5</td></tr><tr><td></td><td>European beech</td><td>352</td><td>17</td><td>91.2 (39.3)</td><td>25.7</td></tr></tbody></table>

Meteorological data were monitored with an automatic meteorological station located in an open field that is 300 m away from the forest site. Soil horizon properties were characterized based on the soil profile description and the measurements carried out by Jonard et al. (2011).

To run the simulations, the values of some model parameters were taken directly from the literature. Other parameters involved in empirical relationships were fitted with either data from previous studies or unpublished monitoring data collected in the study site or in the International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests (ICP Forests) level II plots of Wallonia (Table 2). Potential explanatory variables in Eq. (31) used to estimate height growth were selected by applying a stepwise forward selection procedure based on the Bayesian information criterion (BIC). A multivariate model was then adjusted with the selected variables (Appendix E). The parameters of the npp to gpp ratio relationship, the maintenance respiration per gram of nitrogen at 15 °C, and the PAR use efficiency of sunlit and shaded leaves were adjusted with the nonlinear minimization (nlm) function in R (R Core Team, 2013) based on observed basal area increments (BAIs) using the maximum likelihood approach. This calibration was completed only based on the data of the mixed stand, while the model performances were evaluated with observations from the three stands of the Baileux site.

All the simulations carried out in this study were run with the default option for modelling phenology and water balance (de Wergifosse et al., 2019). In addition, since the tree nutrition and nutrient cycling module was not activated, the tree nutrient status remained constant during the simulations.

## 2.4 Statistical evaluation of model predictions

The quality of the model was evaluated for various combinations of model options (i.e. photosynthesis model of CASTANEA vs. PUE, npp to gpp ratio vs. temperature-dependent maintenance respiration and distance-dependent vs. independent crown extension) by comparing predicted and observed BAIs using several statistical indices and tests such as the normalized average error, the *P* value of the paired *t* test, the regression test, the root mean square error and the Pearson's

correlation (Janssen and Heuberger, 1995). For the regression test, the Deming fitting procedure (mregc function of the mcr package in R) was retained to account for the errors in both the observations and the predictions. For all of the simulations, the water balance module was activated. Some option combinations were therefore not tested, such as the PUE approach without activating the water balance.

The model quality was also evaluated based on its ability to reconstruct the size-growth relationships for the sessile oak and European beech in the three stands in Baileux. The observed and predicted BAIs of the trees (calculated for the 2001–2011 period) were related to their girth at the beginning of the assessment period. A segmented regression was then applied to observations and predictions to determine the girth threshold beyond which the basal area increment (BAI) linearly increases with girth and to estimate the slope of the linear relationship between the BAI and initial girth. The heteroscedasticity of the residuals was accounted for by modelling their standard deviation with a power function of the initial girth. The fitting was carried out using the nlm function in R.

## 2.5 Simulation experiment

To assess how the tree biomass production and its allocation to the different tree compartments were affected by climate conditions and management in the model, we simulated the development of the mixed stand during a dry (2003, with *P* = 948 mm and *T*<sub>air</sub> = 9.88 °C), a normal (2005, with *P* = 1027 mm and *T*<sub>air</sub> = 9.67 °C) and a wet year (2012, with *P* = 1117 mm and *T*<sub>air</sub> = 9.37 °C), and we repeated these simulations after thinning this stand by reducing its basal area by 25%. The biomass production and its allocation were assessed at the stand level as well as at the tree level for seven cohorts (four beech cohorts and three oak cohorts), defined based on the tree species and on the girth-class distribution. For this first simulation experiment, we used the following options: photosynthesis model of CASTANEA, npp to gpp ratio and distance-independent crown extension.

A second simulation experiment was performed to illustrate how the model can be used to predict climate change impacts on the functioning of the forest ecosystem. The growth

<footnote>* Girth at breast height.</footnote>
<table><caption>**Table 2.** Description of model parameters for sessile oak and European beech and origin of their value.</caption><thead><tr><th rowspan="2">Symbol</th><th rowspan="2">Description</th><th rowspan="2">Units</th><th colspan="2">Value</th><th rowspan="2">Origin</th></tr><tr><th>Sessile oak</th><th>European beech</th></tr></thead><tbody><tr><td colspan="6"><strong>Carbon fixation</strong></td></tr><tr><td>k</td><td>extinction coefficient</td><td>m<sup>-1</sup></td><td>0.53</td><td></td><td>fitted with tree growth data of the study site</td></tr><tr><td>PUE<sub>sl</sub></td><td>PAR use efficiency of sunlit leaves</td><td>kg C per mole of photons</td><td>0.00006</td><td>0.000216</td><td>fitted with tree growth data of the study site</td></tr><tr><td>PUE<sub>sh</sub></td><td>PAR use efficiency of shaded leaves</td><td>kg C per mole of photons</td><td>0.00105</td><td>0.000584</td><td>fitted with tree growth data of the study site</td></tr><tr><td colspan="6"><strong>Respiration</strong></td></tr><tr><td>a<sub>sapwood</sub></td><td>parameters of the sapwood area function (a/b/c in Eq. 12)</td><td></td><td>0.00/1.54/0.16</td><td>0.00/0.00/0.52</td><td>fitted with data from André et al. (2010)</td></tr><tr><td>r<sub>npp_gpp</sub></td><td>parameters of the npp to gpp ratio function (α/b in Eq. 8)</td><td></td><td>0.997/-0.386</td><td>0.959/-0.408</td><td>fitted with tree growth data of the study site</td></tr><tr><td>R<sub>Tref</sub></td><td>maintenance respiration per gram of N at the reference temperature (15 °C)</td><td>moles CO<sub>2</sub> per gram of N per hour</td><td>0.000079</td><td>0.000057</td><td>fitted with tree growth data of the study site</td></tr><tr><td>R<sub>gr</sub></td><td>growth respiration per unit biomass increment</td><td>kg C per kg C</td><td>0.2</td><td></td><td>Dufrêne et al. (2005)</td></tr><tr><td>Q<sub>10_leaf/fine root</sub></td><td>temperature dependence coefficient of leaf and fine root respiration</td><td>dimensionless</td><td>2.1</td><td></td><td>Vose and Bolstad (1999)</td></tr><tr><td>Q<sub>10_stem/root</sub></td><td>temperature dependence coefficient of stem and root respiration</td><td>dimensionless</td><td>1.7</td><td></td><td>Epron et al. (2001)</td></tr><tr><td>Q<sub>10ahan</sub></td><td>temperature dependence coefficient of branch respiration</td><td>dimensionless</td><td>2.8</td><td></td><td>Damesin et al. (2002)</td></tr><tr><td colspan="6"><strong>Carbon allocation</strong></td></tr><tr><td>b<sub>leaf</sub></td><td>parameters of the leaf biomass function (α/b/γ in Eq. 15)</td><td>kg C</td><td>0.0026/1.96/1.96</td><td>1.469/2.00/0.00</td><td>Jonard et al. (2006)</td></tr><tr><td>b<sub>structuralAbove</sub></td><td>parameters of the aboveground structural biomass (α/b/γ in Eq. 26)</td><td>kg C</td><td>0.000/263.4/0.969</td><td>0.056/292.8/0.966</td><td>Hounzandj et al. (2015) and Genet et al. (2011)</td></tr><tr><td>r<sub>root_shoot</sub></td><td>root-to-shoot ratio</td><td>kg C per kg C</td><td>0.18</td><td></td><td>Genet et al. (2010)</td></tr><tr><td>r<sub>fr_leaf_min</sub></td><td>minimum fine root to leaf ratio</td><td>kg C per kg C</td><td>0.5</td><td></td><td>literature data compilation</td></tr><tr><td>r<sub>fr_leaf_max</sub></td><td>maximum fine root to leaf ratio</td><td>kg C per kg C</td><td>2.5</td><td></td><td>literature data compilation</td></tr><tr><td>δ<sub>leaf</sub></td><td>leaf relative loss rate</td><td>kg C per kg C per year</td><td>1</td><td></td><td>Grote and Pretzsch (2002)</td></tr><tr><td>δ<sub>fr</sub></td><td>fine root relative loss rate</td><td>kg C per kg C per year</td><td>1</td><td></td><td>Grote and Pretzsch (2002)</td></tr><tr><td>f</td><td>stem form factor</td><td>m<sup>3</sup> per m<sup>3</sup></td><td>0.52</td><td></td><td>Hounzandj et al. (2015) and Genet et al. (2011)</td></tr><tr><td>ρ</td><td>stem volumetric mass</td><td>kg C m<sup>-3</sup></td><td>562.17</td><td>556</td><td>Hounzandj et al. (2015) and Genet et al. (2011)</td></tr><tr><td>rt<sub>leaf</sub></td><td>leaf retranslocation rate</td><td>kg C per kg C per year</td><td>0.4</td><td>0.45</td><td>determined based on tree foliage data of the study site</td></tr><tr><td>r<sub>root</sub></td><td>fine root retranslocation rate</td><td>kg C per kg C per year</td><td>0.4</td><td>0.45</td><td>same values as leaves</td></tr><tr><td>s<sub>branch</sub></td><td>parameters of the branch mortality function (α/b/γ as in Eq. 15)</td><td>kg C</td><td>6.0E-9/3.064/3.064</td><td>5.00E-5/2.681/0.00</td><td>fitted with data from André et al. (2010)</td></tr><tr><td>p<sub>fruit</sub></td><td>parameters of the fruit production function (α/b in Eq. 22)</td><td>kg C</td><td>9.50E-4/2.5</td><td>8.00E-4/2.5</td><td>fitted with litterfall data from level II plots of Wallonia</td></tr><tr><td>dbh<sub>threshold</sub></td><td>threshold dbh for fruit production</td><td>cm</td><td>25</td><td></td><td>field observations</td></tr></tbody></table>
<table><caption>**Table 2. Continued.**</caption><thead><tr><th rowspan="2">Symbol</th><th rowspan="2">Description</th><th rowspan="2">Units</th><th colspan="2">Value</th><th rowspan="2">Origin</th></tr><tr><th>Sessile oak</th><th>European beech</th></tr></thead><tbody><tr><td colspan="6">Tree dimension increment</td></tr><tr><td>hlce%</td><td>height fraction corresponding to the largest crown extension height</td><td>m per m</td><td>0.81</td><td>0.77</td><td>determined based on tree inventory data of the study site</td></tr><tr><td>hcb%</td><td>height fraction corresponding to the crown base height</td><td>m per m</td><td>0.7</td><td>0.61</td><td>determined based on tree inventory data of the study site</td></tr><tr><td>Dd</td><td>parameters of the crown to stem diameter function (α/b/γ/δ in Eq. 10)</td><td>m per m</td><td>16.20/0.0280/0.00/0.00</td><td>10.49/0.00/1379/-2881</td><td>determined based on tree inventory data of the study site</td></tr><tr><td>sh</td><td>coefficient to shift the crown to stem diameter ratio to its maximum</td><td>dimensionless</td><td>1.25</td><td>1.5</td><td>determined based on tree inventory data of the study site</td></tr><tr><td>r<sub>overlapping</sub></td><td>mean crown overlapping ratio</td><td>m per m</td><td>1</td><td>1.2</td><td>determined based on tree inventory data of the study site</td></tr><tr><td>Δhlce<sub>max</sub></td><td>maximum annual change in the largest crown extension height</td><td>m yr<sup>-1</sup></td><td></td><td>0.5</td><td>determined based on tree growth data of the study site</td></tr><tr><td>Δhcb<sub>max</sub></td><td>maximum annual change in the crown base height</td><td>m yr<sup>-1</sup></td><td></td><td>0.5</td><td>determined based on tree growth data of the study site</td></tr></tbody></table>

dynamics in the mixed stand in Baileux was simulated according to three IPCC climate scenarios using the following options: photosynthesis model of CASTANEA, npp to gpp ratio and distance-independent crown extension. The climate scenarios retained for this study were obtained from the global circulation model CNRM-CM5 (Voldoire et al., 2013) based on the representative concentration pathways (RCPs) for atmospheric greenhouse gases described in the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (Collin et al., 2013). The representative concentration pathways (RCP2.6, RCP4.5 and RCP8.5) are characterized by the radiative forcing in the year 2100 relative to preindustrial levels (+2.6, +4.5 and +8.5 W m⁻²). The CNRM-CM5 describes the earth system climate using variables such as air temperature and precipitation on a low-resolution grid (1.4° in latitude and longitude). Although reliable for estimating global warming, such a model fails to capture the local climate variations. Therefore, these climate projections were downscaled by the Royal Meteorological Institute of Belgium (RMI), using the regional climate model ALARO-0 (Giot et al., 2016). The meteorological files that were received from RMI are hourly values of the longwave and shortwave radiation, air temperature, surface temperature, rainfall, specific humidity, zonal and meridional wind speeds and atmospheric pressure with a 4 km spatial resolution. Specific humidity was converted into relative humidity using the Tetens formula (Tetens, 1930). For a reference period (1976–2005), we compared the models predictions with

observed meteorological data and detected some biases, especially for precipitation (overestimation of 27%). To correct these biases, we applied correction factors depending on the month (Maraun and Widmann, 2018). An additive correction factor was used for the bounded variables (radiation, precipitation, relative humidity and wind speed), and a multiplicative was used for the other variables (air and surface temperatures).

For the simulations, two 24-year periods (100 years apart) were considered. The period from 1976 to 1999 served as a historical reference, while the rest of the simulations based on climate projections were conducted for the 2076–2099 period. The simulations were performed by either keeping the CO₂ concentration of the atmosphere constant (i.e. 380 ppm) or allowing it to vary yearly according to the climate scenarios. Each simulation started with the same initial stand (mixed stand in Baileux in 2001) and lasted 24 years; a thinning operation (25 % in basal area) was carried out in 1978 or 2078 and in 1990 or 2090 (12-year cutting cycle). The mean basal area increments obtained with the various climate scenarios were compared using the Tukey multiple comparison test.
<table><caption>**Table 3.** Statistical evaluation of predicted basal area increments (vs. observations) for various combinations of model options using normalized average error (NAE), paired *t* test, regression test, root mean square error (RMSE) or Pearson's correlation (Pearson's *r*). Standard deviation or confidence intervals are provided in parentheses.</caption><thead><tr><th rowspan="2">Model options<br>Tree species</th><th rowspan="2">NAE<br>P value</th><th rowspan="2">Paired t test<br>P value</th><th colspan="2">Orthogonal regression</th><th rowspan="2">RMSE</th><th rowspan="2">Pearson's r</th></tr><tr><th>intercept</th><th>slope</th></tr></thead><tbody><tr><td colspan="7">CASTANEA; npp to gpp ratio; distance-independent crown extension</td></tr><tr><td>European beech</td><td>0.159</td><td>0.00</td><td>1.59(±1.06)</td><td>0.75(±0.05)</td><td>8.64</td><td>0.87</td></tr><tr><td>Sessile oak</td><td>-0.052</td><td>0.18</td><td>5.28(±2.85)</td><td>0.75(±0.14)</td><td>9.33</td><td>0.63</td></tr><tr><td colspan="7">CASTANEA; npp to gpp ratio; distance-dependent crown extension</td></tr><tr><td>European beech</td><td>0.090</td><td>0.04</td><td>1.75(±1.36)</td><td>0.79(±0.07)</td><td>8.95</td><td>0.83</td></tr><tr><td>Sessile oak</td><td>-0.020</td><td>0.61</td><td>4.54(±3.06)</td><td>0.77(±0.14)</td><td>9.11</td><td>0.63</td></tr><tr><td colspan="7">CASTANEA; temperature-dependent maintenance respiration; distance-independent crown extension</td></tr><tr><td>European beech</td><td>0.426</td><td>0.00</td><td>3.55(±2.06)</td><td>0.53(±0.07)</td><td>17.97</td><td>0.74</td></tr><tr><td>Sessile oak</td><td>-0.013</td><td>0.79</td><td>7.18(±3.02)</td><td>0.62(±0.13)</td><td>11.03</td><td>0.59</td></tr><tr><td colspan="7">CASTANEA; temperature-dependent maintenance respiration; distance-dependent crown extension</td></tr><tr><td>European beech</td><td>0.544</td><td>0.00</td><td>2.89(±2.07)</td><td>0.53(±0.06)</td><td>18.60</td><td>0.77</td></tr><tr><td>Sessile oak</td><td>0.054</td><td>0.25</td><td>6.07(±3.42)</td><td>0.64(±0.14)</td><td>11.07</td><td>0.58</td></tr><tr><td colspan="7">PUE; npp to gpp ratio; distance-independent crown extension</td></tr><tr><td>European beech</td><td>0.007</td><td>0.85</td><td>1.60(±1.19)</td><td>0.86(±0.06)</td><td>7.64</td><td>0.85</td></tr><tr><td>Sessile oak</td><td>-0.181</td><td>0.00</td><td>2.93(±4.14)</td><td>1.02(±0.24)</td><td>9.04</td><td>0.54</td></tr><tr><td colspan="7">PUE; npp to gpp ratio; distance-dependent crown extension</td></tr><tr><td>European beech</td><td>-0.110</td><td>0.01</td><td>1.84(±1.57)</td><td>0.96(±0.09)</td><td>8.41</td><td>0.79</td></tr><tr><td>Sessile oak</td><td>-0.223</td><td>0.00</td><td>1.76(±5.61)</td><td>1.16(±0.35)</td><td>9.67</td><td>0.45</td></tr><tr><td colspan="7">PUE; temperature-dependent maintenance respiration; distance-independent crown extension</td></tr><tr><td>European beech</td><td>0.182</td><td>0.01</td><td>4.17(±2.37)</td><td>0.61(±0.09)</td><td>15.43</td><td>0.68</td></tr><tr><td>Sessile oak</td><td>-0.172</td><td>0.00</td><td>4.65(±3.85)</td><td>0.90(±0.21)</td><td>9.31</td><td>0.55</td></tr><tr><td colspan="7">PUE; temperature-dependent maintenance respiration; distance-dependent crown extension</td></tr><tr><td>European beech</td><td>0.223</td><td>0.00</td><td>3.36(±2.47)</td><td>0.64(±0.09)</td><td>14.73</td><td>0.71</td></tr><tr><td>Sessile oak</td><td>-0.176</td><td>0.00</td><td>4.06(±5.02)</td><td>0.95(±0.28)</td><td>9.79</td><td>0.47</td></tr></tbody></table>

# 3 Results

## 3.1 Reconstructed npp vs. predicted gpp

Based on two successive stand inventories (2001 and 2011) and using HETEROFOR in reverse mode (see Sect. 2.2.7), the individual npp was reconstructed and related to the gpp predicted with the photosynthesis method of CASTANEA. The linear relationship between npp and gpp explained 79% and 83% of the variability in sessile oaks and European beeches, respectively (Fig. 2). The intercept was positive and just significantly different from 0 but did not differ between the two trees species. The slope of the relationship was higher for the sessile oak (0.50) than for European beech (0.40).

## 3.2 Model performance in predicting individual basal area increment (BAI)

HETEROFOR was run with different combinations of options for describing photosynthesis (biochemical model of CASTANEA vs. PUE), respiration (npp to gpp ratio vs. temperature-dependent maintenance respiration) and crown extension (distance-dependent vs. distance-independent). The predictions carried out using the photosynthesis routine of CASTANEA were generally slightly better correlated to the observations than those obtained with the PUE approach, which however displayed somewhat lower RMSE (Table 3). For both of the photosynthesis calculation options, the use of the maintenance respiration routine provided less accurate predictions (higher NAE and RMSE and lower Pearson's *r*) than the npp to gpp ratio approach, and the degradation of
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Relationship between the individual npp reconstructed based on successive stand inventories (2001 and 2011) and the gpp predicted with the process-based option (photosynthesis method of CASTANEA) for the three stands. Values in parentheses are 95 % confidence intervals for the intercept and the slope in the equations. The Pearson's correlation between npp and gpp is indicated on the graph.</figcaption></figure>

the model performance due to the maintenance respiration option was more marked for the European beech than for the sessile oak (Table 3). The option for describing crown extension had little effect on prediction quality. Depending on the criterion considered, the options selected for calculating photosynthesis and respiration and the tree species, the distance-independent approach was sometimes the best but not in all cases (Table 3).

For the simulations using the CASTANEA photosynthesis, we retained the npp to gpp approach and the distance-dependent crown extension as the best combination of options since the associated predictions were on average not biased for oak and only slightly biased for beech (Table 3). For this combination of options, the regression of the observed BAIs on the predictions showed however a slight underestimation of the low BAIs and a small overestimation of the high BAIs, which were more pronounced for the European beech than for the sessile oak (Fig. 3). For the PUE method, the npp to gpp ratio and the distance-independent crown extension provided the most accurate predictions (Table 3).

## 3.3 Reconstructing size-growth relationships

The size-growth relationships were very similar between observations and predictions for the mixed stand on which the model was calibrated (Fig. 4). For the European beech in the beech-dominated stand, the predicted increase in BAI with the initial girth was steeper than the observed increase, revealing a slight overestimation of the tree growth (Fig. 4). The proportion of the BAI variance explained by the size-growth relationship ($R^2$) was higher for the European beech than for the sessile oak for both observations and predictions (Fig. 4).

## 3.4 Simulation of climate change impact on tree growth

In the first simulation experiment, the effect of thinning was much more pronounced on the smallest trees than on the biggest ones (Fig. 5). The smallest beech cohort (girth of 0 to 61 cm) almost doubled their annual biomass production after the thinning (+85 %), while the thinning impact on the biggest oak and beech trees was hardly noticeable (+4 % and +2 %, respectively). When looking at the different tree compartments, one may notice that the thinning effect was more pronounced on the structural compartments, i.e. roots, stem and branches (+52 %), than on the functional ones, i.e. fine roots, leaves and fruits (+22 %). While thinning increased the individual biomass production, it decreased the biomass production at the stand level (−15 %).

The biomass production at the stand level was 11 % higher for the normal year than for the dry year (Fig. 5). This effect was observed for all of the cohorts even if it was less marked on the smallest trees (+2 % for the 0–61 cm beech cohort) than on the biggest ones (+13 % for oak and beech trees with a girth larger than 140 cm). Whatever the scale considered (tree or stand), there was nearly no difference in biomass production between the normal and wet year. The climate condition effects were marked only on the structural compartment (+25 %).

When the CO₂ concentration of the atmosphere was fixed, no effect of the climate scenario was detected on stand BAI, but a slight impact was observed on sessile oak BAI, which was higher for the RCP2.6 than for the historical scenario (Fig. 6). For the simulations with a variable atmospheric CO₂ concentration, the differences in total, sessile oak and European beech BAIs were much more pronounced among climate scenarios. For the whole stand as well as for oak and beech separately, BAIs increased the following order (1) historical, (2) RCP2.6, (3) RCP4.5 and (4) RCP8.5, with the stand BAI in each of these RCP scenarios being between 17 % and 72 % higher than that of the historical scenario. All of the scenarios had BAIs that were significantly different from each other, except RCP2.6 and RCP4.5 for the whole stand and the two tree species and historical and RCP2.6 for the European beech (Fig. 6).

## 4 Discussion

Few tree-level, process-based and spatially explicit models have been developed, and these often contain only some of the modules necessary to estimate resource availability (solar radiation, water and nutrients). While descriptions of these models are generally available in the literature, doing an evaluation by making a comparison with tree growth measurements is not always accessible, or such an evaluation has been carried out based on stand-level variables. We have therefore very little information to compare the performances of HETEROFOR at the tree level with those of similar mod-
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of observed and predicted basal area increments (BAIs) for the simulation with the photosynthesis method of CASTANEA, the npp to gpp ratio approach to account for tree respiration and the distance-dependent crown extension (see Table 3). The dashed line represents the Deming regression between observations and predictions, with the shaded area indicating the 95% confidence interval and the solid line the 1:1 relationship.</figcaption></figure>

els. Simioni et al. (2016) faced the same problem with the NOTG 3-D model.

HETEROFOR first estimates the key phenological dates, the radiation interception by trees and the hourly water balance (de Wergifosse et al., 2019). Then, based on the absorbed PAR radiation, individual gpp is calculated with a PUE approach or with the photosynthesis routine of CASTANEA (Dufrêne et al., 2005). Whatever the option retained for calculating tree respiration and crown extension, the photosynthesis routine of CASTANEA and the PUE efficiency approach performed similarly (Table 3). It is quite encouraging that the process-based approach for estimating photosynthesis provided predictions of the same quality as the empirical approach fitted with tree growth data taken on the study site. If no extrapolation to future climate is required, the PUE approach remains however still valuable, especially when hourly meteorological data are lacking. For the three stands in Baileux, we related the npp reconstructed from successive tree inventories with the gpp predicted based on the CASTANEA approach (Fig. 2). The good linear relationships (Pearson's correlation > 0.89) obtained for both oak and beech make us confident in the adaptation of the photosynthesis routine of CASTANEA to the tree level. Furthermore, since the parameters of the photosynthesis routine were taken directly from CASTANEA and not calibrated specifically for HETEROFOR, one can expect that the agreement between the predicted gpp and the reconstructed npp could still be improved.

When comparing the two options available in HETEROFOR for converting gpp into npp, model performances were

systematically better with the npp to gpp ratio approach than with the temperature-dependent routine for the maintenance respiration calculation (Table 3). This can be partly explained because the error in the maintenance respiration calculation results from various sources. At the tree compartment level, uncertainties in the estimation of biomass, sapwood proportion, nitrogen concentration and temperature are multiplied (Eq. 11). Then, the errors made on all of the tree compartments are summed up. Among these uncertainty sources, the inaccuracy in the estimation of the sapwood proportion could explain why the maintenance respiration routine provided better results for the sessile oak than for European beech (Table 3). Since the sapwood of sessile oak can easily be distinguished from the heartwood based on the colour change, we had a lot of sapwood measurements available to fit a relationship. For the European beech, this was not the case; instead, we used a sapwood relationship obtained based on sap flow measurements (Jonard et al., 2011). This relationship could certainly be improved by direct measurements of sapwood made after staining the wood to highlight the living parenchyma. Another way to improve these relationships is to consider the social status of the trees since dominant trees have a higher sapwood depth than the suppressed ones (Rodríguez-Calcerrada et al., 2015). We tried to account for this by estimating the sapwood area based on the tree growth rate, but it did not significantly increase the quality of the predictions. The poor performances obtained with the maintenance respiration option also indicate that the processes at play are still poorly understood and that further research is needed on this topic.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Reconstruction of the size-growth relationships for sessile oak and European beech in the three stands using the photosynthesis method of CASTANEA, the npp to gpp ratio approach to account for tree respiration and the distance-dependent crown extension. The predicted relationships between the individual BAI (calculated for the 2001–2011 period) and the initial girth are compared with observed ones. The solid and dashed lines represent the segmented regression applied, respectively, to observations and predictions to determine the girth threshold beyond which radial growth linearly increases with girth and to estimate the slope of the linear relationship between BAI and initial girth. The 95 % confidence intervals for the intercept and the slope are provided as well as the R² of the model. No relationship was fitted for the European beech in the oak-dominated stand given the lack of data.</figcaption></figure>

The process-based approach for estimating maintenance respiration explicitly accounts for the temperature effect through a $Q_{10}$ function. With the npp to gpp ratio approach, temperature is considered more indirectly by assuming that it affects respiration and photosynthesis in the same proportion, which is valid only in a given range of temperature (< 20 °C) and for non-stressing conditions. Indeed, the optimum temperature for photosynthesis is between 20 and 30 °C, while the optimum temperature for respiration is just below the

temperature of enzyme inactivation (> 45 °C). Therefore, between 30 and 45 °C, photosynthetic rates decrease, but the respiration rate could continue to increase (Yamori et al., 2013). This reasoning however does not consider that the base rate of respiration acclimate to new mean temperature conditions and that this acclimation process tends to maintain the npp to gpp ratio more stably (Collalti and Prentice, 2019). In addition, while water stress reduces both photosynthesis and respiration, its effect on the two processes is not
All cohorts and tree species

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

Sessile oak (circ. 61–105 cm)

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

Sessile oak (circ. 105–140 cm)

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

Sessile oak (circ. 140+ cm)

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

European beech (circ. 0–61 cm)

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

European beech (circ. 61–105 cm)

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

European beech (circ. 105–140 cm)

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

European beech (circ. 140+ cm)

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

**Figure 5.** Effects of climate conditions and thinning on biomass production and its allocation to tree compartment in the mixed stand in Baileux. The data used to make these graphs were obtained by simulations using the following options: photosynthesis model of CASTANEA, npp to gpp ratio and distance-independent crown extension.
Fixed atmospheric CO₂ concentration

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

Variable atmospheric CO₂ concentration

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

**Figure 6.** Basal area increment (BAI) of the mixed stand in Baileux (and of its two main tree species) simulated with climate scenarios produced with the general circulation model CNRM-CM5, downscaled with ALARO-0 and corrected empirically for remaining biases. The simulations were performed by using the CASTANEA method to calculate photosynthesis, the npp to gpp ratio approach and a distance-independent description of crown extension. The CO₂ concentration of the atmosphere was either kept constant (left) or increased with time according to the climate scenario considered (right). Two time periods were considered. The 1976–1999 period was used as a reference for running the model with the historical climate scenario, while the simulations with future climate scenarios were carried out for the 2076–2099 period. The climate scenarios were based on the representative concentration pathways for atmospheric greenhouse gases that are described in the fifth assessment report of IPCC. For a given tree species and CO₂ concentration modality, the scenarios with common letters have BAIs not significantly different from each other (α = 0.05).

necessarily equivalent (Rodríguez-Calcerrada et al., 2014). The main argument in favour of the npp to gpp ratio approach is the tight coupling between respiration and photosynthesis since the substrate for respiration originates from photosynthesis. The npp to gpp ratio is unfortunately neither universal nor constant. It may vary with tree development stage, climate, soil fertility and competition conditions (Collalti and Prentice, 2019). The alternative option based on maintenance respiration calculation is theoretically more appropriate to simulate the impact of climate change, but this is at the expense of less accurate predictions at the tree level. The ideal is to compare the two options to evaluate the prediction uncertainty associated with the modelling of respiration. In the future, the two approaches could be improved. Applying the reconstruction procedure of HETEROFOR on a large diversity of sites would allow us to estimate the npp to gpp ratio in many different situations, to create a function predicting the npp to gpp ratio based on its main drivers and to subsequently use it in the model. In parallel, the respiration calculation could be refined by accounting for thermal acclimation, such as is done in 3D-CMCC (Collalti et al., 2018).

The differences in prediction quality between the two methods of crown extension modelling (distance-dependent approach vs. distance-independent approach) were quite small, probably because the length of the simulation was not sufficient to drastically affect the crown dimensions, which had been initialized based on measurements. Describing mechanisms that govern crown development in interaction with neighbours (mechanical abrasion and crown interpenetration) is however crucial to capturing the nonaddi-

tive effect of species mixtures (Pretzsch, 2014). By accounting for crown plasticity, our distance-dependent approach could help us to better understand how uneven-aged and mixed stands optimize light interception by canopy packing and how they increase productivity (Forrester and Albrecht, 2014; Jucker et al., 2015). To better evaluate the relevance of this approach, the predicted crown development should be compared with precise crown measurements that are repeated over several decades and taken in a large diversity of stand structures. When the model will be calibrated for a larger number of tree species, long-term simulations could also be performed to evaluate to what extent the model is able to reconstruct the empirical relationships that describe tree allometry variations in response to intra- and inter-specific competition. Such relationships were established by del Río et al. (2019) using data from the Spanish National Forest Inventory.

Based on the current evaluation, the process-based variant performs similarly to the more empirical one for photosynthesis and crown extension but not for respiration; this is probably because the processes are better known for photosynthesis. For the best combination of options using the CASTANEA photosynthesis (npp to gpp ratio with distance-dependent crown extension), the Pearson's correlation between measurements and predictions of individual basal increment amounted to 0.83 and 0.63 for European beech and sessile oak, respectively. By comparison, Grote and Pretzsch (2002) obtained a correlation of 0.60 for the individual volume of beech trees with the BALANCE model. This
lower correlation can partly be explained by the integration
of the uncertainty in tree height in the volume estimations.

Individual npp and retranslocated carbon are allocated first to foliage and fine roots and then partitioned between above- and below-ground structural compartments. Based on the derivative and rearrangement of a biomass allometric equation, the increment in aboveground structural biomass is used to estimate the combined increment in dbh and height. This results in a system of one equation with two unknowns (increment in dbh and height). We decided to resolve it by fixing the height growth based on a relationship that takes into account tree size (dbh or height), the height growth potential (height increment if all of the remaining carbon was allocated to height growth) and a light competition index. An intermediate level of sophistication was adopted to describe height growth, between the simple height-dbh allometry and the fine description of tree architecture of functional-structural models. Height-dbh relationships provide a static picture in which age and neighbour effects are confounded and are not suitable to describe individual growth trajectories (Henry and Aarsen, 1999). More sophisticated relationships considering age and dominant height can be used for even-aged stands (Le Moguédec and Dhôte, 2012) but are hardly applicable in uneven-aged stands for which tree age is unknown. On the other hand, the functional-structural models that are based on resource availability at the organ level and use a short time step can only be applied to a limited number of trees given the high computational demand (Letort et al., 2008).

Our individual height growth model was fitted with height
data measured 10 years apart (Appendix E). A large uncer-
tainty was however associated to these data. First, height
measurements were obtained to the nearest metre given the
difficulty to clearly identify the top of the trees in closed
canopy forests. Second, as the height increment was obtained
based on repeated height measurements, the error for this
variable is the sum of those made on the height measure-
ments. Consequently, the uncertainty was more or less of the
same order of magnitude as the expected height growth in
10 years. Despite these uncertainties, a substantial part of the
variability was explained by the model (72% for European
beech; 43% for oak). Among the variables tested, the height
growth potential had the main effect, which is not surpris-
ing since this height growth potential contains the informa-
tion on height increment. We were also able to depict the
effect of light competition. For a same-height growth poten-
tial, trees undergoing stronger light competition seem to in-
vest more carbon for height growth than for dbh increment
(Fig. E1 in Appendix E), which is corroborated by results
of other studies (e.g. Lines et al., 2012). This strategy aims at
minimizing overtopping by neighbours and maximizing light
interception (Jucker et al., 2015). Trouvé et al. (2015) found
similar results and showed the positive effect of stand den-
sity on height growth in the allocation between height and
diameter increment in even-aged stands of sessile oak trees.
The decrease in the red: far red ratio of incident light pro-

motes apical dominance and internode elongation through
the phytochrome system (shade avoidance reaction; Henry
and Aarsen, 1999). By considering the light availability ef-
fect on height growth at the tree level, HETEROFOR adapts
tree allometry to intra- and inter-specific competition, which
is crucial to accounting for mixing effects in structurally
complex stands (del Río et al., 2019).

A first simulation experiment was achieved to assess how
tree biomass production and its allocation to tree compart-
ments respond to climate conditions and thinning. The re-
sults of these simulations are in line with the basic princi-
ples of silviculture; thinning favoured individual tree growth
(especially that of the smaller trees) by redistributing stand
biomass production on a smaller number of trees. At the
stand level, thinning slightly reduced biomass production
since its intensity was substantial and the simulation lasted
only 1 year, which was not sufficient to allow the remain-
ing trees to fill the gaps by extending their crown. Drought
conditions reduced the biomass increment of structural com-
ponents, and this effect was more pronounced on big trees
than on small trees. Indeed, when soil water availability de-
creases, smaller trees maintain a higher stomatal conductiv-
ity because of their lower position in the stand canopy. Func-
tional compartments were less influenced by climate because
carbon is allocated to them in priority in the model. We could
improve the allocation routine by making the fine root to fo-
liage ratio and the root-to-shoot ratio dependent on the mean
soil water availability (Thurm et al., 2017).

We were also quite satisfied to observe that the model
was able to reproduce the size-growth relationship. This ap-
proach describes the growth partitioning among trees in a
stand, which is useful to estimate the mode of competition.
For the three studied stands and the two tree species, the com-
petition was partially size asymmetric, with a resource parti-
tioning in favour of the larger trees (Carl et al., 2018). Within
the studied stands, the European beech trees can be classified
in two groups: a group of small suppressed trees whose ra-
dial growth was close to 0 and which were just surviving and
the rest of the trees (beyond a girth threshold) whose radial
growth linearly increased with girth. Regarding sessile oaks,
nearly all of the trees were in the second group, which can
be related to the fact that sessile oak is a less shade-tolerant
species than European beech. In the this mixed stand, the
nearly perfect match between the predicted and observed re-
lationships indicates that the model was able to reproduce
growth partitioning among trees of different tree species and
size. These very good results can be ascribed to the fact that
the extinction coefficient and the respiration parameters were
calibrated with data from this stand. In the beech-dominated
stand, the model slightly underestimated the radial growth of
the small oak trees and overestimated that of the big beech
trees. In this case, the model seems to allocate too many
resources to the big beech trees which shade the small oak
trees. This could be improved by a model calibration partly
specific to this stand (for the npp to gpp ratio) or by a cali-
bration with data covering a much larger range of stand struc-
tures.

To illustrate one possible application of HETEROFOR, a
second simulation experiment was completed and allowed us
to compare the radial growth predicted for 2076–2099 ac-
cording to three IPCC scenarios with that of the one sim-
ulated for an historical period (1976–1999). When the at-
mospheric CO2 concentration was kept constant (380 ppm),
differences among scenarios remained nonsignificant, except
for sessile oak displaying a slightly higher basal area incre-
ment for the RCP2.6 scenario than for the historical scenario
(Fig. 6). Analysing the model outputs in depth, we found that
this lack of effects resulted from a balance between nega-
tive and positive impacts of climate change. While the in-
crease in air temperature (+0.9 and 3.7 °C for RCP2.6 and
RCP8.5) and in the vegetation period length (+8 and 37 days
for RCP2.6 and RCP8.5) favoured photosynthesis, the more
frequent and intense water stress negatively affected it (de
Wergifosse et al., 2020). The positive and negative effects
of climate change were of the same magnitude for both
tree species and offset each other. For the simulations with
a variable atmospheric CO2 concentration, the differences
among scenarios were much larger, highlighting a strong
CO2 fertilization effect for both sessile oak and European
beech (Fig. 6). These results are in agreement with Reyer et
al. (2014), which used the 4C model to predict productiv-
ity change in Europe according to a large range of climate
change projections. They found npp increases in most Euro-
pean regions (except a few cases in Mediterranean moun-
tains) when considering persistent CO2 effects by using a
variable atmospheric CO2 concentration. Assuming an ac-
climation of photosynthesis to CO2 (by maintaining constant
atmospheric CO2), they predicted increases in northern Eu-
rope, decreases in southern Europe and ambivalent responses
elsewhere in Europe. Similar response patterns were also ob-
tained by Morales et al. (2007). Rötzer et al. (2013) used
the BALANCE model to compare the impact of future and
current climate conditions on the productivity of beeches in
Germany and showed a 30 % decrease in npp without consid-
ering the rise in atmospheric CO2 concentration. After eval-
uating CASTANEA against eddy covariance and tree growth
data in a few highly instrumented sites, Davi et al. (2006)
simulated the trend in gpp and net ecosystem productivity
(NEP) in these sites from 1960 to 2100. For sessile oak and
European beech, they obtained a 53 % and 67 % increase in
gpp and NEP, respectively.

Given the magnitude of the CO₂ fertilization effect (lead-
ing to a 72 % increase in basal area increment in 100 years for
RCP8.5), we conducted retrospective simulations to check
that HETEROFOR reproduces well the increase in produc-
tivity observed by Bontemps et al. (2011) for beech forests
in the north-east of France (data not shown). Based on his-
torical atmospheric CO₂ concentrations, we simulated radial
growth during two periods (1879–1910 vs. 1979–2010) us-
ing the same climate data (obtained by reanalysis for 1979–

2010). These simulations showed a productivity increase of
12% over 100 years. By comparison, Bontemps et al. (2011)
reported productivity increases ranging from 10 to 70% over
100 years, depending on the nitrogen status of the forest. The
increase in radial growth simulated with HETEROFOR for
the mixed stand in Baileux (Fig. 6) seems therefore plaus-
ible but assumes unchanged nutritional status. Increased
productivity generates however higher nutrient demand by
trees, which is not systematically satisfied by larger soil nu-
trient supply, especially in the poorest sites. Consequently,
the augmentation of forest productivity will most likely be
constrained by nutrient availability and give rise to a dete-
rioration of the nutritional status as already observed across
Europe (Jonard et al., 2015). To improve our predictions, nu-
tritional constraints must be taken into account. From this
perspective, a mineral nutrition and nutrient cycle module
was incorporated in HETEROFOR. As it was developed in
parallel to the water balance, some adaptations are needed
for the coupling of the two modules (e.g. change from an
annual to a monthly time step for soil chemistry update). A
complete description and evaluation of the nutrient module
will be provided in a future study.

5 Conclusion and future prospects

Our ambition was to develop a model that is responsive to
both management actions and environmental changes and
would be particularly well adapted to mixed and uneven-aged
stands. We thought that this model had to be tree-level and
spatially explicit and had to consider radiation transfer, water
balance and nutrient cycling with a process-based approach.
Such models were very scarce in the literature. The only ex-
ceptions were BALANCE, iLand and more recently NOTG
3-D. To fill this gap, we elaborated the HETEROFOR model
based on concepts quite different from those used for BAL-
ANCE, iLand and NOTG 3-D. In this study, a first evaluation
of the model performances showed that HETEROFOR was
able to reproduce size-growth relationships in three oak and
beech stands in the Belgian Ardennes. We also noticed that
the npp to gpp ratio option for describing maintenance respi-
ration provides the best results, while the process-based and
empirical approaches perform similarly for photosynthesis
and crown extension. As this model evaluation was limited
to two tree species and one climate, it only provides a first
impression of the model potential.

Here, only the core of HETEROFOR was described. The
water balance and phenology modules are presented and
evaluated in a companion paper (de Wergifosse et al., 2019),
while the nutrient module will be described later. For the
next steps, we plan to couple HETEROFOR with existing
libraries such as the ones for regeneration, genetics and eco-
nomics. As HETEROFOR was developed within the Capsis
platform, it is continually improving thanks to the collabora-
tive dynamics among modellers.
A broader assessment of the model performances will be
carried out based on forest monitoring plots distributed all
over Europe. Indeed, HETEROFOR was designed to be par-
ticularly suitable for the level II plots of ICP Forests. The
processes were described at a scale that facilitates the com-
parison between model predictions and observations. Many
data collected in these plots can be used to initialize and run
the model or to calibrate and evaluate it. HETEROFOR can
also be seen as a tool for integrating forest monitoring data
and quantifying non-measured processes. While it is now cal-
ibrated for oak and beech forests, HETEROFOR will be pa-
rameterized for a large range of tree species in order to use it
for testing and reproducing identity and diversity effects.

Given all of the uncertainties related to climate change im-
pacts, it is unrealistic to believe that a model will accurately
predict the future dynamics of forest growth. However, mod-
els such as HETEROFOR can be very useful to compare sce-
narios. Among others purposes, HETEROFOR can be used
to select the management options that maximize ecosystem
resilience or to quantify uncertainty in the response of forest
ecosystem to climate change.
## Appendix A: Description of the soil heat transfer routine

The temperature of the mineral soil (T; in degrees Celsius) is calculated by soil depth increment (Δz; in metres) using a simplification of the soil heat transfer equation assuming a constant thermal diffusivity (D; in square metres per second) across the soil profile. The thermal diffusivity characterizes the rate of heat transfer within the soil and corresponds to the ratio of the thermal conductivity (K; in watts per metre per kelvin) to the volumetric heat capacity (c_v; in joules per cubic metre per kelvin).

$$ \frac{\partial T}{\partial t} = \frac{1}{c_v} \cdot \frac{\partial}{\partial z} \cdot \left( K \frac{\partial T}{\partial z} \right) \Rightarrow \frac{\partial T}{\partial t} = D \cdot \frac{\partial^2 T}{\partial z^2} \quad (A1) $$

Eq. (A1) can be rewritten according to Anlauf and Liu (1990) and Baker and Don Scott (1998) as follows:

$$ T_{z,t+\Delta t} = T_{z,t} + D \cdot \frac{\Delta t}{\Delta z^2} \cdot (T_{z+\Delta z,t} + T_{z-\Delta z,t} - 2T_{z,t}). \quad (A2) $$

The soil depth increment can be chosen by the user, but it has to be smaller than one-third of the thinnest horizon. The soil depth increment can be slightly modified by the model to ensure that the soil depth is a multiple of the soil depth increment. Then, a stability criterion is checked for each hour, and if it is not respected, the temporal step is divided by two.

$$ K \cdot \frac{\Delta T}{\Delta z^2} < 0.5 \quad (A3) $$

The thermal diffusivity is calculated for each soil horizon based on the thermal conductivity and the volumetric heat capacity and then averaged by weighing according to the horizon thickness. The thermal conductivity is obtained with the empirical model of Kersten (1949) as follows:

$$ K = 0.1442 \cdot (0.9 \cdot \log(\vartheta) - 0.2) \times 10^{0.6243 \rho_b} \quad (\text{for silt or clay soils}) \quad (A4) $$

$$ K = 0.1442 \cdot (0.7 \cdot \log(\vartheta) + 0.4) \times 10^{0.6243 \rho_b} \quad (\text{for sandy soils}), \quad (A5) $$

with ϑ, the gravimetric soil water content (in grams per gram), and ρ_b, the bulk density (in kilograms per cubic metre).

The volumetric heat capacity of soils is approximated through a separation of the soil constituents in solid and liquid phases as follows:

$$ c_v \approx 836 \cdot \rho_b + 4180 \cdot \vartheta \cdot \rho_b \cdot 1000 \cdot \rho_w \quad (A6) $$

with ρ_w, the volumetric mass of water (in kilograms per cubic metre).

To initialize the procedure, the top and bottom temperature during the whole simulation and the initial temperature

at each soil depth must be known. The soil temperature at the top of the mineral soil (just under the forest floor) is given by Eq. (A7) adapted from van Wijk and de Vries (1963) and Chota et al. (2004). The bottom temperature is fixed and corresponds to the mean annual air temperature. This assumption can be made as the soil depth largely exceeds 1 m. The initial temperature is found through a simple interpolation of the temperatures between the soil interface and the bottom.

$$ T_t = \bar{T}_y + \frac{(\bar{T}_{d-1} - \bar{T}_y)}{A_{\text{air}}} \cdot A_{\text{soil}} + \frac{a_{\text{air}}}{2} \cdot \text{red}_d \cdot \sin\left(\omega(t - t_{\text{max}}) + \frac{\pi}{2} - \omega \frac{\Delta z}{\text{Damping}}\right), \quad (A7) $$

with T_y, mean annual air temperature (in degrees Celsius), T_{d-1}, mean air temperature of the previous day (in degrees Celsius), A\_air, annual air temperature amplitude corresponding to the difference between the maximum and the minimum mean daily temperature over the year (in degrees Celsius), A\_soil, a parameter corresponding to the mean annual soil temperature amplitude (in degrees Celsius), a\_air, daily air temperature amplitude (T\_max - T\_min) calculated over the 24 h period centred on the considered time (°C), red\_d, parameter reducing the daily air temperature amplitude to the daily soil temperature amplitude (fixed to 0.13), ω, radial frequency (per hour) is 2π/24, t\_max, hour of the day at which air temperature is maximal (as the sinusoidal shape of the diurnal soil temperature cycle is not perfectly symmetric, t\_max is adapted so that the period between maximum and minimum soil temperature is exactly 12 h), Δz, thickness of organic horizons (in metres), and Damping, a parameter accounting for the phase shift between the diurnal cycle of the air and soil temperature (fixed to 0.0853 after calibration).

The temperature of the organic horizons was obtained as the mean between air temperature and the temperature at the interface between organic horizons and mineral soil.
## Appendix B: Development of Eq. (29)

Equation (29) can be developed in order to isolate $\Delta h$ as follows:

$$ 
\begin{aligned}
\Delta (\mathrm{dh}^2 \cdot h) &= (\mathrm{dh} + \Delta\mathrm{dh})^2 \cdot (h + \Delta h) - \mathrm{dh}^2 \cdot h \\
\Delta (\mathrm{dh}^2 \cdot h) &= (\mathrm{dh}^2 + 2 \cdot \mathrm{dh} \cdot \Delta\mathrm{dh} + (\Delta\mathrm{dh})^2) \cdot (h + \Delta h) - \mathrm{dh}^2 \cdot h \quad (\text{B1})
\end{aligned}
$$

$$ 
\begin{aligned}
\Delta (\mathrm{dh}^2 \cdot h) &= \mathrm{dh}^2 \cdot h + 2 \cdot \mathrm{dh} \cdot \Delta\mathrm{dh} \cdot h + (\Delta\mathrm{dh})^2 \\
& \quad \cdot h + \mathrm{dh}^2 \cdot \Delta h + 2 \cdot \mathrm{dh} \cdot \Delta\mathrm{dh} \\
& \quad \cdot \Delta h + (\Delta\mathrm{dh})^2 \cdot \Delta h - \mathrm{dh}^2 \\
& \quad \cdot h \Delta (\mathrm{dh}^2 \cdot h) \\
&= \Delta\mathrm{dh}^2 \cdot (h + \Delta h) + (\Delta\mathrm{dh})^2 \\
& \quad \cdot (h + \Delta h) + \mathrm{dh}^2 \cdot \Delta h \quad (\text{B2})
\end{aligned}
$$

$$ 
\begin{aligned}
\Delta (\mathrm{dh}^2 \cdot h) &= \Delta h \cdot (\Delta\mathrm{dh}^2 + (\Delta\mathrm{dh})^2 + \mathrm{dh}^2) + h \\
& \quad \cdot (\Delta\mathrm{dh}^2 + (\Delta\mathrm{dh})^2). \quad (\text{B3})
\end{aligned}
$$

Considering $(\Delta\mathrm{dh})^2 \ll \Delta\mathrm{dh}^2 \ll \mathrm{dh}$, an approximation can be done as follows:

$$ \Delta (\mathrm{dh}^2 \cdot h) \approx \Delta h \cdot \mathrm{dh}^2 + h \cdot \Delta\mathrm{dh}^2 \quad (\text{B4}) $$

$$ \Delta h \cdot \mathrm{dh}^2 \approx \Delta (\mathrm{dh}^2 \cdot h) - h \cdot \Delta\mathrm{dh}^2 \quad (\text{B5}) $$

$$ \Delta h \approx \frac{\Delta (\mathrm{dh}^2 \cdot h)}{\mathrm{dh}^2} - \frac{h \cdot \Delta\mathrm{dh}^2}{\mathrm{dh}^2}. \quad (\text{B6}) $$
<table><caption>**Table C1.** Parameters of Eqs. (C3) and (C4) for sessile oak and European beech according to Dagnelie et al. (1999).</caption><thead><tr><th></th><th>a</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th></tr></thead><tbody><tr><td colspan="7">Sessile oak</td></tr><tr><td>C10%</td><td>3.9330</td><td>1.0284</td><td>-0.31611 × 10<sup>-3</sup></td><td>0.44036 × 10<sup>-6</sup></td><td>-0.33113</td><td>-0.28051 × 10<sup>-5</sup></td></tr><tr><td>CR60%</td><td>0.4838</td><td>14.667</td><td>-405.67</td><td></td><td></td><td></td></tr><tr><td colspan="7">European beech</td></tr><tr><td>C10%</td><td>3.8541</td><td>1.0235</td><td>-0.36276 × 10<sup>-3</sup></td><td>0.40063 × 10<sup>-6</sup></td><td>-0.30551</td><td>-0.20411 × 10<sup>-5</sup></td></tr><tr><td>CR60%</td><td>0.5286</td><td>0</td><td>0</td><td></td><td></td><td></td></tr></tbody></table>

**Appendix C: Delevoy height estimation**

The Delevoy height is the height at which stem diameter is half the diameter at breast height and is calculated from taper (in centimetres per metre) as follows:

$$h_{\text{del}} = 1.3 + \frac{\text{dbh} - \frac{\text{dbh}}{2}}{\text{taper}}, \qquad (\text{C1})$$

where the taper is obtained based on the girth at 10% of the tree height (G10%) and the relative girth at 60% of the tree height (RG60%) for which empirical equations are provided by Dagnelie et al. (1999) for several temperate tree species.

$$\text{taper} = \frac{(1 - \text{CR60%}) \cdot \text{C10%}}{0.5 \cdot h \cdot \pi}, \qquad (\text{C2})$$

with

$$\begin{aligned} \text{C10\%} = & a + b \cdot \pi \cdot \text{dbh} + c \cdot (\pi \cdot \text{dbh})^2 + d \cdot (\pi \cdot \text{dbh})^3 \\ & + e \cdot h + f \cdot (\pi \cdot \text{dbh})^2 \cdot h \end{aligned} \qquad (\text{C3})$$

$$\text{CR60\%} = a + \frac{b}{\text{C10\%}} + \frac{c}{\text{C10\%}^2}. \qquad (\text{C4})$$
## Appendix D: Estimation of the height of largest crown extension (hlce) at equilibrium

### Estimation of hlce at equilibrium for a competitor of the same size

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

### Estimation of hlce at equilibrium for a competitor of higher size

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

### Estimation of hlce at equilibrium for a competitor of lower size

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

**Figure D1.** Illustration of the routine used to determine the height of largest crown extension at equilibrium (hlce<sub>eq</sub>) of a target tree in three contrasting situations of competition. A first step consists in determining the intersection between the potential crown of the target tree and the competitor. Then, the hlce<sub>eq</sub> is fixed between the maximum hlce (corresponding to the intersection between potential crowns) and the minimum hlce (which is the height to crown base) based on the relative height of the competitor.
## Appendix E: Height growth modelling results

The main factor explaining the height increment was the so-called height growth potential ($\Delta h_{pot}$), with a quadratic effect for beech and a cubic effect for oak (Table E1 and Fig. E1). For both tree species, the light competition index (LCI) had a negative effect on the height increment, meaning that, for a same-height growth potential, trees under stronger competition for light had a higher height growth than trees within better light conditions. For the European beech, the variable selection procedure selected height (which had a negative effect) to account for tree size, while dbh was retained for sessile oak and had a positive effect. Even if the root-mean-square error was slightly higher for European beech (0.094) than for sessile oak (0.083), the height growth model explained a much larger proportion of the variability in European beeches (72%) than in sessile oaks (43%). This is partly because the height growth range was higher for European beeches.

<table><caption>**Table E1.** Parameters, $R^2$ and RMSE of the height growth model (Eq. 31) for European beech and sessile oak.</caption><thead><tr><th></th><th>European beech</th><th>Sessile oak</th></tr></thead><tbody><tr><td>Intercept</td><td>0.0233</td><td>-0.0562</td></tr><tr><td>dbh (in centimetres)</td><td></td><td>0.0023</td></tr><tr><td>h (in metres)</td><td>-0.0048</td><td></td></tr><tr><td>LCI</td><td>-0.2556</td><td>-0.1874</td></tr><tr><td>(Δdbh<sup>2</sup>h<sub>tot</sub>)/dbh<sup>2</sup> (in metres)</td><td>0.6631</td><td>0.8183</td></tr><tr><td>[(Δdbh<sup>2</sup>h<sub>tot</sub>)/dbh<sup>2</sup>]<sup>2</sup></td><td>-0.1777</td><td>-0.9178</td></tr><tr><td>[(Δdbh<sup>2</sup>h<sub>tot</sub>)/dbh<sup>2</sup>]<sup>3</sup></td><td></td><td>0.4733</td></tr><tr><td>RMSE</td><td>0.09397</td><td>0.083017</td></tr><tr><td>R<sup>2</sup></td><td>0.72</td><td>0.43</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Figure E1.** Effect of the height growth potential on oak and beech height growth for two levels of light competition (strong light competition = light competition index ≤ 0.15, lower light competition = light competition index > 0.15). The solid lines represent the model predictions obtained using Eq. (31) with parameter values of Table E1 and with mean values for dbh, height or the light competition index.</figcaption></figure>
*Code availability.* The source codes of Capsis and HETEROFOR are accessible to all of the members of the Capsis co-development community. Those who want to join this community are welcome but must contact François de Coligny (coligny@cirad.fr) and sign the Capsis charter (http://capsis.cirad.fr/capsis/charter, last access: 29 February 2020). This charter grants access on all the models to the modellers of the Capsis community. The modellers may distribute the Capsis platform with their own model but not with the models of the others without their agreement. Capsis4 is a free software (LGPL licence) which includes the kernel, the generic pilots, the extensions and the libraries. For HETEROFOR, we also choose an LGPL licence and decided to freely distribute it through an installer containing the Capsis4 kernel and the latest version (or any previous one) of HETEROFOR upon request from Mathieu Jonard (mathieu.jonard@uclouvain.be). The version 1.0 used for this paper is available at http://amap-dev.cirad.fr/projects/capsis/files (last access: 29 February 2020, Jonard et al., 2020). The end users can install Capsis from an installer containing only the HETEROFOR model, while the modellers who signed the Capsis charter can access to the complete version of Capsis with all of the models. Depending on your status (end user vs modeller or developer), the instructions to install Capsis are given on the Capsis website (http://capsis.cirad.fr/capsis/documentation, last access: 29 February 2020).

The source code for the modules published in *Geo-scientific Model Development* can be downloaded from https://doi.org/10.5281/zenodo.3591348 (Jonard et al., 2019).

*Data availability.* The data used in this paper are available through the input files for HETEROFOR which are embedded in the installer (see Code availability).

*Supplement.* The supplement related to this article is available online at: https://doi.org/10.5194/gmd-13-905-2020-supplement.

*Author contributions.* MJ, FA, FdC, LdW, NB and HD developed the model code. FA carried out the calibration. MJ, FA and LdW performed the simulations and analysed the model outputs. MJ prepared the paper with contributions from all coauthors.

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* We are grateful to the RMI (Royal Meteorological Institute of Belgium) for providing us climate projection data for the Baileux site. We also would like to thank the four anonymous reviewers for their constructive comments which helped us to significantly improve the quality of the paper.

*Financial support.* This research has been supported by the Service Public de Wallonie (SPW/DGO 3/DNF) (Accord-Cadre de Recherche et Vulgarisation Forestières 2014–2019), the Fonds de

la Recherche Scientifique – FNRS (PDR-WISD grant no. 09; the SustainFor project) and the Fonds de la Recherche Scientifique – FNRS, Fonds pour la Formation à la Recherche dans l'Industrie et dans l'Agriculture (grant no. 1.E005.18).

*Review statement.* This paper was edited by Min-Hui Lo and reviewed by Rüdiger Grote and three anonymous referees.

## References

* Aber, J., Neilson, R. P., McNulty, S., Lenihan, J. M., Bachelet, D., and Drapek, R. J.: Forest processes and global environmental change: predicting the effects of individual and multiple stressors, Bioscience, 51, 735–751, 2001.

* André, F., Jonard, M., and Ponette, Q.: Precipitation water storage capacity in a temperate mixed oak–beech canopy, Hydrol. Process., 22, 4130–4141, 2008a.

* André, F., Jonard, M., and Ponette, Q.: Influence of species and rain event characteristics on stemflow volume in a temperate mixed oak–beech stand, Hydrol. Process., 22, 4455–4466, 2008b.

* André, F., Jonard, M., and Ponette, Q.: Biomass and nutrient content of sessile oak (*Quercus petraea* (Matt.) Liebl.) and beech (*Fagus sylvatica* L.) stem and branches in a mixed stand in southern Belgium, Sci. Total Environ., 408, 2285–2294, 2010.

* André, F., Jonard, M., Jonard, F., and Ponette, Q.: Spatial and temporal patterns of throughfall volume in a deciduous mixed-species stand, J. Hydrol., 400, 244–255, 2011.

* Anlauf, R. and Liu, Y. P.: Simulation of simple one-dimensional transport processes, in: Models for Processes in the Soil: Programs and Exercises, edited by: Richter, J., Catena Verlag, Cremlingen-Desdett, Germany, 1990.

* Baker, D. L. and Don Scott, H.: A limited tutorial on using finite differences in soil physics problems, available at: http://www.aquarien.com/findif/Findifa4.html (last access: 29 February 2020), 1998.

* Ball, J. T., Woodrow, I. E., and Berry, J. A.: A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions, in: Progress in Photosynthesis Research, edited by: Biggins, J., vol. 4, 221–224, 1987.

* Berger, U., Piou, C., Schiffers, K., and Grimm, V.: Competition among plants: concepts, individual-based modelling approaches, and a proposal for a future research strategy, Perspect. Plant Ecol., 9, 121–135, 2008.

* Boisvenue, C. and Running, S. W.: Impacts of climate change on natural forest productivity–evidence since the middle of the 20th century, Glob. Change Biol., 12, 862–882, 2006.

* Bontemps, J. D., Hervé, J. C., Leban, J. M., and Dhôte, J. F.: Nitrogen footprint in a long-term observation of forest growth over the twentieth century, Trees, 25, 237–251, 2011.

* Bosc, A.: EMILION, a tree functional-structural model: presentation and first application to the analysis of branch carbon balance, Ann. Forest Sci., 57, 555–569, 2000.

* Bravo, F., Fabrika, M., Ammer, C., Barreiro, S., Bielak, K., Coll, L., Fonseca, T., Kangur, A., Löf, M., Merganicova, K., Pach, M., Pretzsch, H., Stojanovic, D., Schuler, L., Peric, S., Rötzer, T., del Rio, M., Dodan, M., and Bravo-Oviedo,
A.: Modelling approaches for mixed forests dynamics prognosis, Research gaps and opportunities, Forest Syst., 28, eR002, doi:10.5424/fs/2019281-14342, 2019.

Campioli, M., Vincke, C., Jonard, M., Kint, V., Demarée, G., and Ponette, Q.: Current status and predicted impact of climate change on forest production and biogeochemistry in the temperate oceanic European zone: review and prospects for Belgium as a case study, J. Forest Res., 17, 1–18, 2012.

Cantarello, E., Newton, A. C., Martin, P. A., Evans, P. M., Gosal, A., and Lucash, M. S.: Quantifying resilience of multiple ecosystem services and biodiversity in a temperate forest landscape, Ecol. Evol., 7, 9661–9675, 2017.

Carl, C., Biber, P., Veste, M., Landgraf, D., and Pretzsch, H.: Key drivers of competition and growth partitioning among Robinia pseudoacacia L. trees, Forest Ecol. Manage., 430, 86–93, 2018.

Charlton, S. R. and Parkhurst, D. L.: Modules based on the geochemical model PHREEQC for use in scripting and programming languages, Comput. Geosci., 37, 1653–1663, 2011.

Cichota, R., Elias, E. A., and van Lier, Q. D. J.: Testing a finite-difference model for soil heat transfer by comparing numerical and analytical solutions, Environ. Model. Softw., 19, 495–506, 2004.

Collalti, A., Trotta, C., Keenan, T. F., Ibrom, A., Bond-Lamberty, B., Grote, R., Vicca, S., Reyer, C. P. O., Migliavacca, M., Veroustraete, F., Anav, A., Campioli, M., Scoccimarro, E., Šigut, L., Grieco, E., Cescatti, A., and Matteucci G.: Thinning can reduce losses in carbon use efficiency and carbon stocks in managed forests under warmer climate, J. Adv. Model. Earth Sy., 10, 2427–2452, 2018.

Collalti, A. and Prentice, I. C.: Is NPP proportional to GPP? Warning's hypothesis 20 years on, Tree Physiol. 39, 1473–1483, 2019.

Collins, M., Knutti, R., Arblaster, J., Dufresne, J. L., Fichefet, T., Friedlingstein, P., Gao, X., Gutowski, W. J., Johns, T., Krinner, G., Shongwe, M., Tebaldi, C., Weaver, A. J., and Wehner, M.: Long-term Climate Change: Projections, Commitments and Irreversibility, in: Climate Change 2013: The Physical Science Basis, Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G. K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P. M., Cambridge University Press, Cambridge, UK, New York, NY, USA, 2013.

Courbaud, B., de Coligny, F., and Cordonnier, T.: Simulating radiation distribution in a heterogeneous Norway spruce forest on a slope, Agr. Forest Meteorol., 116, 1–18, 2003.

Dagnelie, P., Palm, R., Rondeux, J., and Thill, A.: Tables de cubage des arbres et des peuplements forestiers, Les Presses Agronomiques de Gembloux, Gembloux, Belgique, 1999.

Damesin, C., Ceschia, E., Le Goff, N., Ottorini, J. M., and Dufrêne, E.: Stem and branch respiration of beech: from tree measurements to estimations at the stand level, New Phytol., 153, 159–172, 2002.

Davi, H., Dufrêne, E., Francois, C., Le Maire, G., Loustau, D., Bosc, A., Rambal, S., Granier, A., and Moors, E.: Sensitivity of water and carbon fluxes to climate changes from 1960 to 2100 in European forest ecosystems, Agr. Forest Meteorol., 141, 35–56, 2006.

Davi, H., Cailleret, M., Restoux, G., Amm, A., Pichot, C., and Fady, B.: Disentangling the factors driving tree reproduction, Ecosphere, 7, e01389, https://doi.org/10.1002/ecs2.1389, 2016.

del Río, M., Bravo-Oviedo, A., Ruiz-Peinado, R., and Condés, S.: Tree allometry variation in response to intra- and inter-specific competitions, Trees, 33, 121–138, 2019.

de Wergifosse, L., André, F., Beudez, N., de Coligny, F., Goosse, H., Jonard, F., Ponette, Q., Titeux, H., Vincke, C., and Jonard, M.: HETEROFOR 1.0: a spatially explicit model for exploring the response of structurally complex forests to uncertain future conditions. II. Phenology and water cycle, Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2019-201, in review, 2019.

de Wergifosse, L., André, F., Goosse, H., Caluwaerts, S., De Cruz, L., De Troch, R., and Van Schaeybroeck, B.: CO₂ fertilization, water deficit and vegetation period drive the response of mixed broadleaved forests to a changing climate in Wallonia, Ann. Forest Sci., in review, 2020.

Dufour-Kowalski, S., Courbaud, B., Dreyfus, P., Meredieu, C., and de Coligny, F.: Capsis: an open software framework and community for forest growth modelling, Ann. Forest Sci., 69, 221–233, 2012.

Dufrêne, E., Davi, H., François, C., Le Maire, G., Le Dantec, V., and Granier, A.: Modelling carbon and water cycles in a beech forest: Part I: Model description and uncertainty analysis on modelled NEE, Ecol. Model., 185, 407–436, 2005.

Duursma, R. A. and Medlyn, B. E.: MAESPA: a model to study interactions between water limitation, environmental drivers and vegetation function at tree and stand levels, with an example application to [CO₂] × drought interactions, Geosci. Model Dev., 5, 919–940, https://doi.org/10.5194/gmd-5-919-2012, 2012.

Ennos, R., Cottrell, J., Hall, J., and O'Brien, D.: Is the introduction of novel exotic forest tree species a rational response to rapid environmental change? – A British perspective, Forest Ecol. Manag., 432, 718–728, 2019.

Epron, D., Le Dantec, V., Dufrêne, E., and Granier, A.: Seasonal dynamics of soil carbon efflux and simulated rhizosphere respiration in a beech forest, Tree Physiol., 21, 145–152, 2001.

Farquhar, G. D., von Caemmerer, S., and Berry, J. A.: A biochemical model of photosynthetic CO₂ assimilation in leaves of C3 species, Planta, 149, 78–80, 1980.

Fernandez-Martinez, M., Vicca, S., Janssens, I. A., Sardans, J., Luyssaert, S., Campioli, M., Chapin III, F. S., Ciais, P., Malhi, Y., Obersteiner, M., Papale, D., Piao, S. L., Reichstein, M., Rodà, F., and Peñuelas, J.: Nutrient availability as the key regulator of global forest carbon balance. Nat. Clim. Change, 4, 471–476, 2014.

Forrester, D. I. and Albrecht, A. T.: Light absorption and light-use efficiency in mixtures of Abies alba and Picea abies along a productivity gradient, Forest Ecol. Manag., 328, 94–102, 2014.

Genet, H., Bréda, N., and Dufrêne, E.: Age-related variation in carbon allocation at tree and stand scales in beech (*Fagus sylvatica* L.) and sessile oak (*Quercus petraea* (Matt.) Liebl.) using a chronosequence approach, Tree Physiol., 30, 177–192, 2010.

Genet, A., Wernsdörfer, H., Jonard, M., Pretzsch, H., Rauch, M., Ponette, Q., Nys, C., Legout, A., Ranger, J., Vallet P., and Saint-André, L.: Ontogeny partly explains the apparent heterogeneity of published biomass equations for Fagus sylvatica in central Europe, Forest Ecol. Manag., 261, 1188–1202, 2011.

Giot, O., Termonia, P., Degrauwe, D., De Troch, R., Caluwaerts, S., Smet, G., Berckmans, J., Deckmyn, A., De Cruz, L., De Meutter, P., Duerinckx, A., Gerard, L., Hamdi, R., Van den Bergh, J., Van Ginderachter, M., and Van Schaeybroeck, B.: Validation
of the ALARO-0 model within the EURO-CORDEX framework, Geosci. Model Dev., 9, 1143–1152, https://doi.org/10.5194/gmd-9-1143-2016, 2016.

Greene, D., Messier, C., Asselin, H., and Fortin, M.: The effect of light availability and basal area on cone production in Abies bal- samea and Picea glauca, Can. J. Bot., 80, 370–377, 2002.

Grote, R. and Pretzsch, H.: A model for individual tree development based on physiological processes, Plant Biol., 4, 167–180, 2002.

Haxeltine, A. and Prentice, I. C.: A general model for the light use efficiency of primary production by terrestrial ecosystems, Funct. Ecol., 10, 551–561, 1996.

Helmisaari, H. S., Derome, J., Nöjd, P., and Kukkola, M.: Fine root biomass in relation to site and stand characteristics in Nor- way spruce and Scots pine stands, Tree Physiol., 27, 1493–1504, 2007.

Henry, H. A. L. and Aarssen, L. W.: The interpretation of stem diameter-height allometry in trees: biomechanical constraints, neighbour effects, or biased regressions?, Ecol. Lett., 2, 89–97, 1999.

Hounzandji, A. P., Jonard, M., Nys, C., Saint-André, L., and Ponette, Q.: Improving the robustness of biomass functions: from empirical to functional approaches, Ann. Forest Sci., 72, 795–810, 2015.

Janssen, P. H. M. and Heuberger, P. S. C.: Calibration of process-oriented models, Ecol. Model., 83, 55–66, 1995.

Jonard, M., André, F., and Ponette, Q.: Modeling leaf dispersal in mixed hardwood forests using a ballistic approach, Ecology, 87, 2306–2318, 2006.

Jonard, M., André, F., Jonard, F., Mouton, N., Procès, P., and Ponette, Q.: Soil carbon dioxide efflux in pure and mixed stands of oak and beech, Ann. Forest Sci., 64, 141–150, 2007.

Jonard, M., André, F., and Ponette, Q.: Tree species mediated effects on leaf litter dynamics in pure and mixed stands of oak and beech, Can. J. Forest Res., 38, 528–538, 2008.

Jonard, F., André, F., Ponette, Q., Vincke, C., and Jonard, M.: Sap flux density and stomatal conductance of European beech and sessile oak trees in pure and mixed stands during the summer drought of 2003, J. Hydrol., 409, 371–381, 2011.

Jonard, M., Fürst, A., Verstraeten, A., Thimonier, A., Timmermann, V., Potocic, N., Waldner, P., Benham, S., Hansen, K., Merila, P., Ponette, Q., De La Cruz, A. C., Roskams, P., Nicolas, M., Croisé, L., Ingerslev, M., Matteuci, G., Decinti, B., Bascietto, M., and Rautio, P.: Tree mineral nutrition is deteriorating in Europe, Global Change Biol., 21, 418–430, 2015.

Jonard, M., André, F., and de Wergifosse, L.: Code of HETEROFOR 1.0, https://doi.org/10.5281/zenodo.3591348, 2019.

Jonard, M., André, F., and de Wergifosse, L.: Installer of HETEROFOR 1.0, http://amap-dev.cirad.fr/projects/capsis/files, last access: 29 February 2020.

Jucker, T., Bouriaud, O., and Coomes, D. A.: Crown plasticity enables trees to optimize canopy packing in mixed-species forest, Funct. Ecol., 29, 1078–1086, 2015.

Kersten, M. S.: Laboratory research for the determination of the thermal properties of soils, University of Minnesota, Engineering Experiment Station, 1949.

Le Moguédec, G. and Dhôte, J. F.: Fagacées: a tree-centered growth and yield model for sessile oak (*Quercus petraea* L.) and common beech (*Fagus sylvatica* L.), Ann. Forest Sci., 69, 257–269, 2012.

Letort, V., Cournède, P. H., Mathieu, A., de Reffye, P., and Constant, T.: Parametric identification of a functional-structural tree growth model and application to beech trees (*Fagus sylvatica*), Funct. Plant Biol., 35, 951–963, 2008.

Ligot, G., Balandier, P., Courbaud, B., Jonard, M., and Kneeshaw, D.: Managing understory light to maintain a mixture of species with different shade tolerance, Forest Ecol. Manag., 327, 189–200, 2014.

Lindner, M., Maroschek, M., Netherer, S., Kremer, A., Barbati, A., Garcia-Gonzalo, J., Seidl, R., Delzon, S., Corona, P., Kolström, M., lexer, M. J., and Marchetti, M.: Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems, Forest Ecol. Manag., 259, 698–709, 2010.

Lindner, M., Fitzgerald, J. B., Zimmermann, N. E., Reyer, C., Delzon, S., van der Maaten, E., Schelhaas, M.-J. Lasch, P., Eggers, J., van der Maaten-Theunissen, M., Suckow, F., Psomas, A., Poulter, B., and Hanewinkel, M.: Climate change and European forests: What do we know, what are the uncertainties, and what are the implications for forest management? J. Environ. Manage., 146, 69–83, 2014.

Lines, E. R., Zavala, M. A., Purves, D. W., and Coomes, D. A.: Predictable changes in aboveground allometry of trees along gradients of temperature, aridity and competition, Glob. Ecol. Biogeogr., 21, 1017–1028, 2012.

Mäkelä, A.: Implications of the pipe model theory on dry matter partitioning and height growth in trees, J. Theor. Biol., 123, 103–120, 1986.

Mäkelä, A. and Valentine, H. T.: The ratio of NPP to GPP: evidence of change over the course of stand development, Tree Physiol., 21, 1015–1030, 2001.

Manion, P. D.: Tree Disease Concepts, Prentice-Hall, Englewood Cliffs, NJ, 402 pp., 1981.

Maraun, D. and Widmann, M.: Statistical downscaling and bias correction for climate research, Cambridge University Press, Cambridge, 2018.

Mellert, K. H. and Göttlein, A.: Comparison of new foliar nutrient thresholds derived from van den Burg's literature compilation with established central European references, Eur. J. Forest Res., 131, 1461–1472, 2012.

Messier, C., Puettmann, K., Chazdon, R., Andersson, K. P., Angers, V. A., Brotons, L., Filotas, E., Tittler, R., Parrrott, L., and Levin, S. A.: From management to stewardship: viewing forests as complex adaptive systems in an uncertain world, Conserv. Lett., 8, 368–377, 2015.

Monteith, J. L.: Climate and the efficiency of crop production in Britain, Philos. Trans. Royal Soc. B, 281, 277–294, 1977.

Morales, P., Hickler, T., Rowell, D. P., Smith, B., and Sykes, M. T.: Changes in European ecosystem productivity and carbon balance driven by regional climate model output, Global Change Biol., 13, 108–122, 2007.

Oliver, T. H., Heard, M. S., Isaac, N. J. B., Roy, D. B., Procter, D., Eigenbrod, F., Freckleton, R., Hector, A., Orme, C. D. L., Petchey, O. L., Proença, V., Raffaelli, D., Suttle, K. B., Mace, G. M., Martín López, B., Woodcock, B. A., and Bullock, J. M.: Biodiversity and resilience of ecosystem functions, Trends Ecol Evol., 30, 673–684, 2015.

Penning de Vries, F. W. T.: The cost of maintenance processes in plant cells, Ann. Bot., 39, 77–92, 1975.
Pretzsch, H.: Canopy space filling and tree crown morphology in mixed-species stands compared with monocultures, Forest Ecol. Manag., 327, 251–264, 2014.

Pretzsch, H., Forrester, D. I., and Rötzer, T.: Representation of species mixing in forest growth models, A review and perspective, Ecol. Model., 313, 276–292, 2015.

R Core Team: R: A language and environment for statistical computing, 2013.

Reyer, C., Lasch-Born, P., Suckow, F., Gutsch, M., Murawski, A., and Pilz, T.: Projections of regional changes in forest net primary productivity for different tree species in Europe driven by climate change and carbon dioxide, Ann. Forest Sci., 71, 211–225, 2014.

Rodríguez-Calcerrada, J., Martin-St. Paul, N. K., Lempereur, M., Ourcival, J. M., del Carmen del Rey, M., Joffre, R., and Rambal, S.: Stem CO₂ efflux and its contribution to ecosystem CO₂ efflux decrease with drought in a Mediterranean forest stand, Agr. Forest Meteorol., 195–196, 61–72, 2014.

Rodríguez-Calcerrada, J., López, R., Salomón, R., Gordaliza, G. G., Valbuena-Carabaña, M., Oleksyn, J., and Gil, L.: Stem CO₂ efflux in six co-occurring tree species: underlying factors and ecological implications, Plant Cell Environ., 38, 1104–1115, 2015.

Rötzer, T., Liao, Y., Goergen, K., Schüler, G., and Pretzsch, H.: Modelling the impact of climate change on the productivity and water-use efficiency of a central European beech forest, Clim. Res., 58, 81–95, 2013.

Ryan, M. G.: Effects of climate change on plant respiration, Ecol. Appl., 1, 157–167, 1991.

Ryan, M. G. and Yoder, B. J.: Hydraulic limits to tree height and tree growth. What keeps trees from growing beyond a certain height? BioScience, 47, 235–242, 1997.

Schäfer, K. V. R., Oren, R., and Tenhunen, J. D.: The effect of tree height on crown level stomatal conductance, Plant Cell Environ., 23, 365–375, 2000.

Schwalm, C. R. and Ek, A. R.: A process-based model of forest ecosystems driven by meteorology, Ecol. Model., 179, 317–348, 2004.

Simioni, G., Marie, G., and Huc, R.: Influence of vegetation spatial structure on growth and water fluxes of a mixed forest: Results from the NOTG 3D model, Ecol. Model., 328, 119–135, 2016.

Teh, C. B. S.: Introduction to mathematical modeling of crop growth. How the equations are derived and assembled into a computer program, Brown Walker Press, Boca Raton, Florida, USA, 2006.

Tetens, O.: Uber einige meteorologische Begriffe, Z. Geophys., 6, 297–309, 1930.

Thompson, I., Mackey, B., McNulty, S., and Mosseler, A.: Forest resilience, biodiversity, and climate change. A synthesis of the biodiversity/resilience/stability relationship in forest ecosystems. Secretariat of the Convention on Biological Diversity, Montreal, Technical Series no. 43, 2009.

Thurm, E. A., Biber, P., and Pretzsch, H.: Stem growth is favored at expenses of root growth in mixed stands and humid conditions for Douglas-fir (Pseudotsuga menziesii) and European beech (Fagus sylvatica), Trees-Struct. Funct., 31, 349–365, 2017.

Trouvé, R., Bontemps, J. D., Seynave, I., Collet, C., and Lebourgeois, F.: Stand density, tree social status and water stress influence allocation in height and diameter growth of *Quercus petraea* (Liebl.), Tree Physiol., 35, 1035–1046, 2015.

Van Wijk, W. R. and De Vries, D. A.: Periodic temperature variations in a homogenous soil, in: Physics of Plant Environment, edited by: Van Wijk, W. R., North Holland Publishing Company, Amsterdam, 1963.

Villar, R., Held, A. A., and Merino, J.: Dark leaf respiration in light and darkness of an evergreen and a deciduous plant species, Plant Physiol., 107, 421–427, 1995.

Voldoire, A., Sanchez-Gomez, E., Salas y Mélia, D., Decharme, B., Cassou, C., Sénési, S., Valcke, S., Beau, I., Alias, A., Chevallier, M., Déqué, M., Deshayes, J., Douville, H., Fernandez, E., Madec, G., Maisonnave, E., Moine, M. P., Planton, S., Saint-Martin, D., Szopa, S., Tyteca, S., Alkama, R., Belamari, S., Braun, A., Coquart, L., and Chauvin, F.: The CNRM-CM5.1 global climate model: description and basic evaluation, Clim. Dynam., 40, 2091–2121, 2013.

Vose, J. M. and Bolstad, P. V.: Challenges to modelling NPP in diverse eastern deciduous forests: species-level comparisons of foliar respiration responses to temperature and nitrogen, Ecol. Model., 122, 165–174, 1999.

Yamori, W., Hikosaka, K., and Way, D. A.: Temperature response of photosynthesis in C₃, C₄, and CAM plants: temperature acclimation and temperature adaptation, Photosynth. Res., 119, 101–117, 2013.