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

Abstract


Introduction
Forest structure and composition result from soil and climate conditions, management and natural disturbances.All 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 make no doubt, their magnitude and the way they will occur locally remain largely uncertain (Lindner et al., 2014).Designing silvicultural systems and selecting tree species adapted to future conditions seems therefore a risky bet (Ennos et al., 2019).Messier et al. (2015) proposed another approach recognizing that forests are complex adaptive systems whose future dynamics is inherently uncertain.To maintain the ability of forests to provide a large range of goods and services whatever the future conditions, their resilience and adaptability must be improved by favouring 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 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 would provide results only in the long run given the life span of trees.Scenario analysis based on model simulations are therefore useful to select the most promising management strategies and to evaluate their long-term sustainability.To explore forest response to new silvicultural practices and to unexperienced climate conditions in a realistic way, one needs new process-based models able to deal with mixed and structurally complex stands and to incorporate uncertainties in future conditions (Pretzsch et al., 2015).
In connection with the traditional forestry viewing forests as a stable system 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 eco-physiological models to better understand short-and long-term forest ecosystem response 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 eco-physiological 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 for 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 9 at the stand, 11 at the cohort and 16 at the tree level.While cohort models allow to describe 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 position.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 eco-physiological experiments (Duursma and Medlyn, 2012).MAESPA is however not suitable for multi-year simulations since it contains no routine for carbon allocation, respiration and tree dimensional growth.EMILION is also restricted to one-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 eco-physiological concepts compared to MAESPA and photosynthesis is calculated with a 10-day 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.The approach used for modelling nutrient cycling is however very simple.Soil is not partitioned into horizons and the soil chemistry processes (e.g.ion exchange, mineral weathering) are not described although they are essential to estimate bioavailability of the major nutrients other than N (P, K, Mg, Ca).Not considered in the review of Pretzsch et al. (2015), iLand is another individual-based model describing the eco-physiological processes with an intermediate level of detail using simplified eco-physiological concepts (such as the radiation use efficiency approach) in order simulate forest dynamics also at the landscape scale.Later, Simioni et al. (2016) developed the NOTG 3D 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 (a few years) rather 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 FORrests.While the BALANCE and iLand model 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 based on other concepts.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 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 water dynamics, rainfall partitioning when passing through the canopy), our model rests on a simpler representation of tree structure than BALANCE.Second, we aimed at incorporating a detailed tree nutrition and nutrient cycling module since we realized the necessity to integrate nutritional constraints 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 in a collaborative modelling platform dedicated to tree growth and stand dynamics.Among the various platforms, CAPSIS was the only one allowing multi-model integration and providing a user-friendly interface (Dufour-Kowalski et al., 2012).HETEROFOR was therefore progressively elaborated through the integration of various modules (light interception, phenology, water cycling, photosynthesis and respiration, carbon allocation, mineral nutrition and nutrient cycling) within CAPSIS.The advantage of such a platform is to use common development environment, model execution system, userinterface and visualization tools and to share data structures, objects, methods and libraries.5 To simulate the response of forests to management and changing environmental conditions, integrate and structure the existing knowledge into process-based models is essential but not sufficient.These models must also be documented and evaluated in order to know exactly their strengths and limits when analysing their outputs.The objectives of this paper are (i) to describe the carbon-related processes of HETEROFOR (photosynthesis, respiration, carbon allocation and tree dimensional growth), (ii) evaluate the model ability in reconstructing tree growth in three broadleaved stands of contrasted species composition and 10 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 various IPCC climate scenarios.

Overall operation of the HETEROFOR model
HETEROFOR is a model integrated in the CAPSIS (Computer-Aided Projection for Strategies in Silviculture) platform dedicated to forest growth and dynamics modelling (Dufour-Kowalski et al., 2012).CAPSIS provides to HETEROFOR the execution system and the methods necessary to run simulations and display the results.When running simulations with HETEROFOR, CAPSIS creates a new project in which the variables describing the forest state are stored at a yearly time step, starting from the initial forest characteristics (initial step).Though some data structures and methods are shared with other models integrated in CAPSIS, the initialisation and evolution procedures are specific to HETEROFOR.
For the initialization, HETEROFOR loads a series of files containing tree species parameters, input data on tree (location, 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.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 successive yearly time step, the phenological periods for each deciduous species (leaf development, leaf colouring and shedding) are defined for the next step from meteorological data.When no meteorological measurements are available, the vegetation period is defined by the user who provides the budburst and the leaf shedding dates.Knowing the key phenological dates and the rates of leaf expansion, colouring and falling, the foliage state of the deciduous species is predicted at any time during the year (de Wergifosse et al., in review).It is characterized by the proportions of leaf biomass and of green leaves relatively to complete leaf development, which are key variables to simulate energy, water and carbon fluxes within the forest ecosystem.The proportion of green leaves impacts photosynthesis, leaf respiration and tree transpiration, as these processes are not active anymore on discoloured leaves which however still intercept solar radiation and rainfall.Based on a ray tracing approach, the SAMSARALIGHT library of CAPSIS (Courbaud et al., 2003) calculates the solar radiation absorbed by the trunk and the crown of each individual tree and the radiation transmitted to the ground.This allows HETEROFOR to estimate the proportions of incident radiation absorbed by the trunk and the crown of each tree and the part transmitted to the ground either on average over the whole vegetation period (simplified budget) or hourly for several key dates (detailed budget).These proportions and the incident radiation measured in the meteorological station are used during the next step to compute the hourly global, direct and diffuse radiation absorbed per unit bark or leaf area.Predicting how solar energy is distributed within the forest ecosystem is necessary to estimate foliage, bark and soil evaporation, tree transpiration and leaf photosynthesis.
Every hour, HETEROFOR performs a water balance and updates the water content of each horizon.Rainfall is partitioned in throughfall, stemflow and interception (Andre et al., 2008a;2008b and2011).Part of the rainfall reaches directly the ground (throughfall) while the rest is intercepted by foliage and bark.These two tree compartments both have a certain water storage capacity which is regenerated by evaporation.When the foliage is saturated, the overflow joins the throughfall flux whose proportion increases.As the bark saturates, water flows along the trunk to form stemflow.Foliage and bark storage capacity as well as stemflow proportion are determined at the tree level and then upscaled to the stand level, while evaporation from these surfaces is evaluated at the stand scale.Throughfall is also determined at the stand level as the difference between incident rainfall and the abovementioned fluxes.Throughfall and stemflow supply the first soil horizon (forest floor) with water while soil evaporation and root uptake deplete it.The water evaporation from the soil (as well as from the foliage and the bark) is calculated at stand scale with the Penman-Monteith equation based on the solar radiation absorbed by each component.Using the same equation, individual tree transpiration is estimated by determining the stomatal conductance from tree characteristics, soil extractable water and meteorological conditions.The distribution of root water uptake among the soil horizons is done according to the water potential and the vertical distribution of fine roots.Water exchanges between soil horizons are considered as water inputs (capillary rise) or outputs (drainage).This soil water transfers are calculated based on the water potential gradients according to the Darcy law and using pedotransfer functions to determined soil hydraulic properties.All these soil water fluxes are considered at the stand level (de Wergifosse et al., in review).
The gross primary production of each tree (gpp) is either obtained based on a radiation use efficiency approach distinguishing sunlit and shaded leaves or calculated hourly using the Farquhar et al. (1980) model.The latter is analytically coupled to the stomatal conductance model proposed by Ball et al. (1987).The photosynthesis is computed using the Library CASTANEA 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 radiation (PAR) absorbed per unit leaf area and the mean soil water potential.gpp is then converted to net primary production (npp) after subtraction of growth and maintenance respiration.Maintenance respiration is either considered as a proportion of gpp (depending on the crown to stem diameter ratio) or calculated for each tree compartment by considering the living biomass, the nitrogen concentration and a Q10 function for the temperature dependency following Ryan (1991) as in Dufrêne et al. (2005).Carbon allocation is made in priority to foliage and fine roots by ensuring a functional balance between carbon fixation and nutrient uptake through a fine root to leaf biomass ratio depending on the tree nutritional status (Helmisaari et al., 2007).Allometric relationships are then used to describe carbon allocation to structural components (trunk, branches and structural roots) and to derive tree dimensional growth (diameter at breast height, total height, height to crown base, height of largest crown extension, crown radii in 4 directions) while considering competition with neighbouring trees (Fig. 1).
Knowing the chemical composition of the tree compartments for a given tree nutrient status, HETEROFOR computes the individual tree nutrient requirements based on the estimated growth rate and deduces the tree nutrient demand after subtraction of the amount of re-translocated nutrient.On another hand, the potential nutrient uptake 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 maintenance respiration and through the effect of the tree nutrient status on fine root allocation.
The soil chemistry is characterized at the 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.This compartment receives the nutrients coming from atmospheric deposition, organic matter mineralization and primary mineral weathering, and is depleted by root uptake and immobilization in micro-organisms.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 processes regulating the carbon fluxes (Fig. 1) while the phenology and water balance modules are presented in a companion paper (de Wergifosse et al., in review) and the nutrient cycling and tree nutrition module will be described later in a third paper.

Initialization
To initialize HETEROFOR, the relative position (x, y, z) and the main dimensions of each tree must be provided: girth at breast height (gbh in cm), height (h in m), height of maximum crown extension (hlce in m), height to crown base (hcb in m) and crown radii in the four cardinal directions (cr in m).During the initialization phase, the biomass of each tree compartment is calculated according to the equations 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: a file describing soil horizon properties and another one for the hourly meteorology.

Gross primary production
The annual gross primary production of each tree (gpp in kgC yr -1 ) 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).Whatever the option retained, a series of 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 MJ yr - 1 ).aRAD is then converted in PAR (aPAR in mol photons yr -1 ) by supposing that 46% of the solar radiation (RAD) is PAR and that 1 MJ is equivalent to 4.55 moles of photons.The diffuse and direct components of aPAR are also considered (aPARdiff and aPARdir in mol photons yr -1 ).While all the leaves receive diffuse PAR, only sunlit leaves absorb direct PAR.To estimate the sunlit leaf proportion (Propsl) at the tree level, HETEROFOR uses an adaptation of the classical stand-scale approach based on the Beer-Lambert law (Teh, 2006) with , the extinction coefficient, , the leaf area index (m² m -2 ).
At the individual scale, the leaf area index is calculated by dividing the tree leaf area (aleaf in m 2 ) by the crown projection area (cpa in m²).The value obtained is then multiplied by the light competition index (LCI in MJ MJ -1 ) to account for the shading effect of the neighbouring trees: where LCI is the ratio between the absorbed radiation calculated with and without neighbouring trees in SAMSARALIGHT.LCI ranges from 1 (no light competition) to 0 (no light reaching the tree).
To adapt the PAR use efficiency concept (PUE) at the tree level, we considered a distinct PUE for sunlit (sl) and shaded (sh) leaves and calculated an average PUE weighted as follows: This pue is then used to calculate gpp based on aPAR and a reducer accounting for water stress ( #2 3$45) ): The default value of #2 3$45) is 1 but, when the hydrological module is activated, it is set to the ratio between the actual and the potential (i.e., considering no soil water limitation) tree transpiration (9 $:4;$ and 9 +*4 , in l per year).This ratio estimates the fraction of the vegetation period during which stomata are partially or totally closed due to 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.#2 3$45) = 4 <= 4 >< • ?
(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 with the soil water availability characterized in HETEROFOR based on the mean soil water potential (see de Wergifosse et al., in review).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).
The photosynthesis method 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 (µmol photons m -2 s -1 ) by the proportion of direct PAR absorbed by sunlit leaves.For a tree, this proportion is by default fixed 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 mol photons.m -².yr -1 ) and the incident PAR cumulated over the same period (in mol photons m -² yr -1 ).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 phenological stage.

Growth and maintenance respiration
gpp is converted to annual net primary production (npp in kgC yr -1 ) 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).@ = 6 • A++_C++ D2 (6) @ = 6 − F − 6 (7) Mäkelä and Valentine (2001) showed that the npp to gpp ratio changes with tree size.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 size (height, dbh, crown radius, crown volume, crown to stem diameter ratio, aboveground volume or biomass) on the npp to gpp ratio.The best relationship was obtained with the crown to stem diameter ratio (Dd in m m -1 ) 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.
As the crown to stem diameter ratio changes during the course of the tree development for some tree species, we standardized it to obtain a crown to stem diameter index (D2 @2#G).
where H and Jare parameters and D2 @2#G is defined as : with D2, the crown to stem diameter ratio determined from the tree mean crown radius (? O5$A in m) and diameter at breast height (dbh in m), D2 +)5N , the crown to stem diameter ratio predicted based on the girth at breast height (gbh in cm): In Eq. ( 7), maintenance respiration is calculated for each tree by summing the maintenance respiration of each organ estimated from the nitrogen content of its living biomass and considering a Q10 function for the temperature dependency.During daytime, the inhibition of foliage respiration by light is taken into account by considering that this inhibition reduces respiration by 62% (Villard et al., 1995).
with P *)C$A , the organ biomass (kg of organic matter), Y Z[ZAC , the fraction of living biomass, \]^, the nitrogen concentration (g kg -1 ), 8 _ 1 , the maintenance respiration per g of N at the reference temperature (15°C), T, is the air temperature for aboveground organs 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 organs.The sapwood proportion is derived from the sapwood area (7 $+3**N in cm²) determined based on an empirical function of organ diameter (∅ *)C$A in cm): Growth respiration is the sum of the organ growth respiration which is proportional to the organ biomass increment (see sect. 2.2.4): where 8 C) is the growth respiration per unit biomass increment (kgC kgC -1 ).

Carbon allocation and dimensional growth
For each tree, the npp and the carbon retranslocated from leaves and roots ( 9 5$l and 9 lZA5 )**4 in kgC yr -1 ) are distributed among the various tree compartments at the end of the year.9 5$l and 9 lZA5 )**4 are determined as follows : where P 5$l and P lZA5 )**4 are the tree leaf and fine root biomasses (kgC), U 5$l and U lZA5 )**4 are the leaf and fine root turnover rates (kgC kgC -1 yr -1 ), and 9 5$l and 9 lZA5 )**4 are the leaf and fine root retranslocation rates (kgC kgC -1 ).P 5$l is estimated with an allometric equation based on the stem diameter at breast height (dbh in cm) and on the crown to stem diameter ratio (Dd): 15) P lZA5 )**4 is deduced from the leaf biomass using the fine root to leaf ratio ( lZA5 )**4 4* l* Z$C5 ): lZA5 )**4_ 5$l takes a value between a minimum ( lZA5 )**4_ 5$l_OZA ) and maximum ( lZA5 )**4_ 5$l_O$o ) ratio depending on the tree nutritional status, in accordance with the concept of functional balance (Mäkela 1986).This means that a higher ratio is used (more carbon allocation to fine roots) 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) > 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 on the optimum and deficiency thresholds (Mellert and Göttlein, 2012).p979"q ]"9 r#@9 = \s* Z$) t;4)Z5A4^ M5lZ:Z5A:u v+4ZO;O M5lZ:Z5A:u The leaf and fine root litter amounts (q 5$l and q lZA5 )**4 in kgC yr -1 ) are estimated based on the turnover rate taking into account the retranslocation: In the allocation, priority is given to leaves and fine roots.The carbon allocated to leaves corresponds to the annual leaf production ( 5$l in kgC yr -1 ) which is equal to the amount of leaves fallen the previous year plus the leaf biomass change (∆P 5$l in kgC yr -1 ): where ∆P 5$l is determined by : ∆P 5$l = P 5$l 4 − P 5$l 4 (20) with P 5$l 4 and P 5$l 4 being 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: where P lZA5 )**4 4 is provided by Eq. ( 16).
When the carbon allocated to leaf and fine root 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 ( l);Z4 in kgC yr -1 ) is estimated with an allometric equation similar to Eq. ( 15) and is considered directly proportional to the light competition index.A threshold dbh (2Pℎ 4T)5 T* N in cm) is fixed below which no fruit production occurs.
Part of the carbon is also used to compensate for branch and root mortality.The branch mortality (q S)$A:T in kgC yr -1 ) is described with an equation of the same form as Eq. ( 15) while the structural root mortality (q )**4 in kgC yr -1 ) 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 organ growth: ∆P 4);:4;)$ = @ + 9 − 5$l − lZA5 )**4 − l);Z4 − q S)$A:T − q )**4 (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 ( )**4_ T**4 ): ∆P 4);:4;)$ _S5 *3 = ∆P 4);:4;)$ − ∆P 4);:4;)$ _$S*[5 The increment in aboveground structural biomass is then used to determine the combined increment in dbh and total height (h in m) based on an allometric equation used to predict aboveground woody biomass (Genet et al., 2011;Hounzandji et al., 2015): Deriving this equation and rearranging terms gives: The development of the left term provides: which can be further developed (see Appendix B for details) to isolate ∆ℎ: From Eq. ( 30), we know that the height increment can be expressed as a function of ∆i NST V •Tk NST V .In the following, we refer to it as the height growth potential (∆ℎ +*4 ) since it corresponds to the height increment if all the remaining carbon was allocated to height growth.Contrary to the other term of Eq. ( 30) NST V } which is unknown, this height growth potential can be evaluated at this step by dividing the result of Eq. ( 28) by dbh².However, depending on the level of competition for light and on the tree size, only part of this height growth potential will be effectively realised for 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): The dbh increment is then determined by rearranging Eq. ( 29): The increments in root, stem and branch biomass are obtained as follows: with f is the form coefficient (m 3 m -3 ), " is the stem volumetric mass (kgC m -3 ), ℎ N5 is the Delevoy height (m) corresponding to the height at which stem diameter is half the diameter at breast height The branch and root biomasses are then distributed in 3 categories defined based on the diameter: small branches/roots < 4 cm, medium branches/roots between 4 and 7 cm, coarse branches/roots > 7 cm.The proportions of small, medium and coarse branches/roots are determined based on the equations of Hounzandji et al. (2015).

Crown extension
Depending on whether the competition with the neighbouring trees is taken into account or not, the crown dynamics can be describe by two different approaches.When local competition is not considered (distance-independent approach), change in crown dimensions are derived from dbh or height increment based on empirical relationships: ∆ℎ"?# = ℎ"?#% • ∆ℎ (36) ∆ℎ?P = ℎ?P% • ∆ℎ (37) where ℎ?P% and ℎ"?#% are the proportions of the total height corresponding to the height to crown base (ℎ?P in m) and to the height of largest crown extension (ℎ"?# in m), respectively; ∆? is the change in crown radius (in m) whatever the direction; D2 +)5N 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 4 sectors according to the 4 cardinal directions (North between 315° and 45°,East between 45° and 135°,South between 135° and 225°,West between 225° and 315°).In each sector, the tree which 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., two times the maximal crown radius: 10 m), no competitor is considered.For each main direction, the model calculates an hlce at equilibrium (ℎ"?# 5-in m) for the target tree.This hlce at equilibrium is located between a minimum (ℎ?P in m) and a maximum (ℎ"?# O$o in m).ℎ"?# O$o 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 (? +*4 in m, see below) and to the crown length (ℎ − ℎ?P).ℎ"?# 5-is positioned between the minimum and the maximum values according to the competition intensity estimated based on the target tree and the competitor heights (ℎ 4$)C54 and ℎ :*O+ in m) as well as the hcb of the target tree (Appendix D): ℎ"?# 5-= ℎ?P + ℎ"?# O$o − ℎ?P • F7G X0, Fr@ 1, The four values of ℎ"?# 5-are then averaged (ℎ"?# 5-_O5$A ).
where ∆ℎ"?# O$o is the maximum change in ℎ"?# allowed by the model.
The change in the four crown radii is calculated based on crown radii at equilibrium (? 5-in m) which are estimated by considering the competitive strength of the target and neighbouring trees.For a given direction, ?5-is calculated based on the potential (free growth) crown radius of the target tree (? +*4_4$)C54 in m) and of its competitor (? +*4_:*O+ in m), the distance between the two trees (d in m) and the crown overlap ratio ( *[5) $+ in m m -1 ): The potential crown radius (? +*4 ) of a tree if determined by: where Ddpred is the crown to stem diameter ratio estimated by Eq. ( 10) and sh is a coefficient allowing to shift from the mean to the maximum Ddpred.
The crown overlap ratio is estimated by considering neighbouring trees of the same species two by two and by 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.

Tree harvesting and mortality
During the simulation, thinning can be achieved at each annual step either (i) by selecting the trees from a list or a map or according to tree characteristics (tree species, age, dbh, height,…), or (ii) by defining the number of trees to be thinned per diameter class using an interactive histogram, or (iii) by 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 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: where P 5$l and P 5$l_:*)) 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).

Growth reconstruction
HETEROFOR was adapted to allow the user to run it in reverse mode starting from the known increments in dbh and h to reconstruct 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 requires 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 use the model equations to reconstruct each step and evaluate among other individual npp.The npp is obtained by re-arranging Eq. ( 23) in which the carbon allocated to the structural biomass is calculated from the dbh and h increments using Eq. ( 27), ( 25) and ( 24).The carbon allocated to leaf, fine root and fruit production is determined respectively with Eq. ( 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).

Input variables and parameter setting for a case study
The model was tested in three stands contrasting in structure and species composition.These stands were located close to each other (< 1km) 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 oak (Quercus petraea Liebl.) and European beech (Fagus sylvatica L.) and lies on acid brown earth soil (luvisol according to the FAO soil taxonomy) with a moder humus and an AhBwC profile.The soil has been developed on a loamy and stony solifluxion sheet in which weathering products of the bedrock (Lower Devonian: sandstone and schist) were mixed with added periglacial loess.
By the end of the 19 th 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 beech.
In 2001, the area was covered by even-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(Jonard et al., , 2007(Jonard et al., , 2008;;André et al., 2008aAndré et al., , 2008bAndré et al., , 2010André et al., , 2011)): two plots were located in stands dominated either by sessile oak or by 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, crown diameters in two directions) at the end of the years 2001 and 2011.
Meteorological data were monitored with an automatic meteorological station located in an open field 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 either with data from previous studies or with unpublished monitoring data collected in the study site or in the ICP Forests level II plots of Wallonia (Table 2).Potential explanatory variables of 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 g of N at 15°C and the PAR use efficiency of sunlit and shaded leaves were adjusted with the nlm function of R (R Core Team, 2013) based on observed basal area increments (BAIs) using the maximum likelihood approach.This fitting was achieved only based on the data of the mixed stand while the model performances were evaluated on the tree stands.

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, 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 (Janssens and Heuberger, 1995).For the regression test, the Deming fitting procedure (mcreg function of the mcr package in R) was retained to account for the errors on both the observations and the predictions.
The model quality was also evaluated based on its ability to reconstruct the size -growth relationships for sessile oak and European beech in the three stands of 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 BAI linearly increases with girth and to estimate the slope of the linear relationship between BAI and initial girth.The heteroscedasticity of the residuals was accounted by modelling their standard deviation with a power function of the initial girth.The fitting was carried out using the nlm function in R.

Simulation experiment
To illustrate how the model can be used to predict climate change impacts on forest ecosystem functioning, the growth dynamics in the mixed stand of Baileux was simulated according to three IPCC climate scenarios using the following options: The CNRM-CM5 describes the earth system climate using variables such as air temperature and precipitations on a lowresolution 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 radiations, 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 precipitations (overestimation of 27%).The biases were corrected by adding to the predictions (or by multiplying them with) a correction factor specific to the month (Maraun and Widmann, 2018).An additive correction factor was used for the bounded variables (radiations, precipitation, relative humidity, wind speed) and a multiplicative one 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 either by keeping the CO2 concentration of the atmosphere constant (i.e., 380 ppm) or by allowing it to vary according to the years and climate scenarios.Each simulation started with the same initial stand (mixed stand of Baileux in 2001) and lasted 24 years; a thinning operation (25% in basal area) was achieved in 1978 or 2078 and in 1990 or 2090 (12-year cutting cycle).The mean basal area increment obtained with the various climate scenarios were compared using the Tukey multiple comparison test.

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 for sessile oak and for European beech, 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 sessile oak (0.50) than for European beech (0.40).

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 -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 options of photosynthesis calculation, the use of the maintenance respiration routine provided less accurate predictions (higher NAE and RMSE and lower Person's r) than the npp to gpp ratio approach and the degradation of the model performance due to the maintenance respiration option was more marked for European beech than for Sessile oak (Table 3).The option for describing crown extension had little effect on prediction quality.Depending on the criterion considered, on the options selected for calculating photosynthesis and respiration and on the tree species, the distance-independent approach was sometimes the best but not systematically (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 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 European beech than for Sessile oak (Fig. 3).For the PUE method, the npp to gpp ratio and the distanceindependent crown extension provided the most accurate predictions (Table 3).

Sessile oak European beech
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.

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 one revealing a slight overestimation of the tree growth (Fig. 4).The proportion of the BAI variance explained by the size -growth relationship (R²) was higher for European beech than for sessile oak for both observations and predictions (Fig. 4).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.

Simulation of climate change impact on tree growth
When the CO2 concentration of the atmosphere was fixed, no effect of the climate scenario was detected on stand BAI but a slight impact was assessed on sessile oak BAI which was higher for the RCP2.6 than for the historical scenario (Fig. 5).For the simulations with a variable atmospheric CO2 concentration, the difference in total, sessile oak and European beech BAI were much more pronounced between climate scenarios.For the whole stand as well as for oak and beech separately, BAI increased in the order -historical, RCP2.6,RCP4.5 and RCP8.5 -, with the stand BAI of these RCP scenarios being between 17 and 72% higher than that of the historical scenario.All scenarios had a BAI 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 European beech (Fig. 5).

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 a description of these models is generally available in the literature, their evaluation by comparison with tree growth measurements is not always accessible or was carried out based on stand-level variables.We have therefore very few information to compare the performances of HETEROFOR at the tree level with those of similar models.Simioni et al. (2016) faced the same problem with the NOTG 3D model.
HETEROFOR first estimates the key phenological dates, the radiation interception by trees and the hourly water balance (de Wergifosse et al., in review).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).This is quite encouraging that the process-based approach for estimating photosynthesis provided predictions of the same quality than 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 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 for 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 maintenance respiration calculation (Table 3).This can be explained since 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 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 sessile oak than for European beech (Table 3).Since the sapwood of sessile oak can easily be distinguish from the heartwood based on the colour change, we had a lot of sapwood measurements available to fit a relationship.For 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 one (Rodriguez-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 increased the quality of the predictions.The poor performances obtained with the maintenance respiration option also indicates that the processes at play are still poorly understood and that further research are needed on this topic.
The process-based approach for estimating maintenance respiration accounts explicitly for the temperature effect through a Q10 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 respiration rate could continue to increase (Yamori et al., 2013).In addition, while water stress reduces both photosynthesis and respiration, its effect on the two processes is not necessarily equivalent (Rodriguez-Calcerrada et al., 2014).Compared to the npp to gpp ratio approach, the maintenance respiration calculation seems 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.
Alternatively, one can choose one or the other option depending on the aim of the simulation.In the future, we could also implement a third option accounting for thermal acclimation such as in 3D-CMCC (Collalti et al., 2018).
The differences in prediction quality between the two methods of crown extension modelling (distance-dependent vsindependent 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.In addition, describing mechanisms that governs crown development in interaction with neighbours (mechanical abrasion, crown interpenetration) is indeed crucial to capture non-additive effect of species mixtures (Pretzsch, 2014).By accounting for crown plasticity, our distance-dependent approach could help better understand how uneven-aged and mixed stands optimize light interception by canopy packing and how they increase productivity (Forester and Albrecht, 2014;Juncker et al., 2015).To fully evaluate the interest of this approach, the predicted crown development should be compared with precise crown measurements repeated over several decades and taken in a large diversity of stand structures.
Based on the current evaluation, the process-based variant perform similarly than the more empirical one for photosynthesis and crown extension but not for respiration, probably because the processes are better known for photosynthesis.For the best combination of options using the CASTANEA photosynthesis (npp to gpp ratio/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 on tree height in the volume estimations.The HETEROFOR performances in terms of tree growth are quite good and could still be improved by a Bayesian calibration of the parameters.
Individual npp and retranslocated C are allocated first to foliage and fine roots and then partitioned between above-and belowground 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 taking into account tree size (dbh or height), the height growth potential (height increment if all 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 functionalstructural models.Height-dbh relationships provides 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 based on resource availability at organ level and using a short time step can only be applied to a limited number of trees given the long computing times (Letort et al., 2008).
Our individual height growth model was fitted with height data measured ten years apart (Appendix E).A large uncertainty was however associated to these data.First, height measurements were obtained to the nearest meter given the difficulty to clearly identify the top of the trees in closed canopy forests.Second, two errors were added since the height increment was obtained by the difference between two height measurements.Consequently, the uncertainty was more or less of the same order of magnitude than the expected height growth in ten 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 surprising since this height growth potential noisily contains the information on height increment.We were also able to highlight the effect of light competition.For a same height growth potential, trees undergoing stronger light competition seem to invest 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 density on height growth in the allocation between height and diameter increment in even-aged stands of sessile oak.The decrease in the red:far red ratio of incident light promotes apical dominance and internode elongation through the phytochrome system (shade avoidance reaction, Henry and Aarsen, 1999).
We were also quite satisfied to observe that the model was able to reproduce the size-growth relationship.This relationship is characterized by two parameters: the threshold which defines the girth beyond which radial growth linearly increases with girth and the slope providing the theoretical maximum growth efficiency (Le Moguédec and Dhôte, 2012).For European beech, the observed threshold was 49 cm and was easy to detect visually since there were many trees with a girth inferior to that exhibiting nearly no basal area increment.For sessile oak, the observed threshold was more variable than for European beech (41 and 60 cm) and nearly all trees had a higher girth.This can be related to the fact that sessile oak is a less shadetolerant species than European beech.
To illustrate one possible application of HETEROFOR, a simulation experiment was achieved and allowed us to compare the radial growth predicted for 2076-2099 according to three IPCC scenarios with that simulated for an historical period .When atmospheric CO2 concentration was kept constant (380 ppm), differences among scenarios remained non-significant, except for sessile oak displaying a slightly higher basal area increment for the RCP2.6 than for the historical scenario (Fig. 5).Analyzing in-depth the model outputs, we found that this lack of effects resulted from a balance between negative and positive impacts of climate change.While the increase in air temperature (+0.9 and 3.7°C for RCP2.6 and 8.5) and in the vegetation period length (+8 and 37 days for RCP2.6 and 8.5) favoured photosynthesis, the more frequent and intense water stress negatively affected it (de Wergifosse et al., in review).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. 5).These results are in agreement with Reyer et al. (2014) who used the 4C model to predict productivity change in Europe according to a large range of climate change projections.They found NPP increases in most European regions (except a few cases in Mediterranean mountains) when considering persistent CO2 effects by using variable atmospheric CO2 concentration.Assuming an acclimation of photosynthesis to CO2 (by maintaining atmospheric CO2 constant), they predicted increases in Northern, decreases in Southern and ambivalent responses elsewhere in Europe.Similar response patterns were also obtained 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 beech in Germany and showed a 30% decrease in NPP without considering the rise in atmospheric CO2 concentration.After evaluating CASTANEA against eddy covariance and tree growth data in a few highly instrumented sites, Davi et al. (2006) simulated the trend in GGP 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 CO2 fertilization effect (leading 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 productivity observed by Bontemps et al. (2011) for beech forests in the north-east of France (data not shown).Based on historical atmospheric CO2 concentrations, we simulated radial growth during two periods (1879-1910 vs 1979-2010) using the same climate data (obtained by re-analysis 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. 5) seems therefore plausible but assumes unchanged nutritional status.Increased productivity generates however higher nutrient demand by trees, which is not systematically satisfied by larger soil nutrient 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 deterioration of the nutritional status as already observed across Europe (Jonard et al., 2015).To improve our predictions, nutritional constraints must be taken into account.In 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 a perfect 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.

Conclusion and future prospects
Our ambition was to develop a model responsive to both management actions and environmental changes that would be particularly well adapted to mixed and uneven-aged stands.We thought that this model had to be tree-level, spatially explicit and process-based.Except BALANCE, iLand and more recently NOTG 3D, no such a model existed in the literature and took simultaneously the radiation transfer, the water balance and the nutrient budget into account.To fill this gap, we elaborated the HETEROFOR model based on concepts quite different from those used for BALANCE, iLAND and NOTG 3D.In this study, a first evaluation of the model performances showed that HETEROFOR predicts well individual radial growth and is able to reproduce size-growth relationships.We also noticed that the most empirical option for describing maintenance respiration provides the best results while the process-based and empirical approaches perform similarly for photosynthesis and crown extension.
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., in review) while the nutrient module will be described later.For the next steps, we plan to couple HETEROFOR with existing libraries such as regeneration, genetics and economics.As HETEROFOR was developed within the CAPSIS platform, it is continually improving thanks to the collaborative 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 particularly suitable for the level II plots of ICP Forests.The processes were described at a scale that facilitates the comparison 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 calibrated for oak and beech forests, HETEROFOR will be parameterised for a large range of tree species in order to use it for testing and reproducing identity and diversity effects.
Given all the uncertainties related to climate change impacts, it is an illusion to believe that a model will predict accurately the future dynamics of forest growth.However, models such HETEROFOR can be very useful to compare scenarios.Among others, HETEROFOR can be used to select the management options that maximise ecosystem resilience or to quantify uncertainty in the response of forest ecosystem to climate change.

Code availability
The source code of CAPSIS and HETEROFOR is accessible to all 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)or Nicolas Beudez (nicolas.beudez@inra.fr)and sign the CAPSIS charter (http://capsis.cirad.fr/capsis/charter).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 license 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. 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 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).

Data availability
The data used in this paper are available through the input files for HETEROFOR which are embedded in the installer (see sect. 6).
with T º » , mean annual air temperature (°C), T º L , mean air temperature of the previous day (°C), A ÇÅÏ , annual air temperature amplitude corresponding to the difference between the maximum and the minimum mean daily temperature over the year (°C), A ÃÄÅAE , parameter corresponding to the mean annual soil temperature amplitude (°C), a ÇÅÏ , daily air temperature amplitude T ÌÇ − T ÌÅÍ calculated over the 24 hour period centered on the considered time (°C), red L , parmeter reducing the daily air temperature amplitude to the daily soil temperature amplitude (fixed to 0.13) ω, radial frequency (h -1 ) = hË h-, t • É¿Ê , 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 • É¿Ê is adapted so that the period between maximum and minimum soil temperature is exactly 12 hours), ∆z, thickness of organic horizons (m), Damping, 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 C -Delevoy height estimation
The Delevoy height is the height at which stem diameter is half the diameter at breast height and is calculated as follows from taper (cm m -1 ): 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:

Appendix E -Height growth modelling results
The main factor explaining the height increment was the so-called height growth potential (∆ℎ +*4 ) with a quadratic effect for beech and a cubic effect for oak (Table E1, Fig. E1).For both tree species, the light competition index (LCI had a negative effect on 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 European beech, the variable selection procedure led to 5 select 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 for European beech (72%) than for sessile oak (43%), partly because the height growth range was higher for European beech.

Figure 1 .
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 or with a biochemical model of photosynthesis (coupling with the CASTANEA library).The photosynthesis calculation depends on the soil water content which is updated hourly thanks to the water balance module described in details in de Wergifosse et al. (in review).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 these processes (radiation interception, photosynthesis and respiration as well as 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 organs (roots, trunk and branches) using a fixed proportion for the below-ground part.dbh and height growth (∆dbh, ∆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 -independent approach.
OZA and ∆? O$o being respectively the minimum and the maximum change in ?allowed by the model.They are obtained similarly as ?+*4 : ∆cr O photosynthesis model, 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 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 RCP8.5) are characterized by the radiative forcing in the year 2100 relative to preindustrial levels (+2.6 W m -2 , +4.5 W m -2 , +8.5 W m -2 ).

Figure 2 .
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.

Figure 4 .
Figure 4. Reconstruction of the size -growth relationships for sessile oak and European beech in the three stands using the photosynthesis method of CASTANEA 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

Figure 5 .
Figure5.Basal area increment (BAI) of the mixed stand in Baileux (and of its two main tree species) simulated with climate scenarios produced with the GCM 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 CO2 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.1976-1999 was used as a reference period for running the model with the historical climate scenario while the simulations with future climate scenarios were achieved for the 2076-2099 period.The climate scenarios were based on the representative concentration pathways for atmospheric greenhouse gases described in the fifth assessment report of IPPC.For a given tree species and CO2 concentration modality, the scenarios with common letters have a BAI not significantly different from each other (α=0.05).

Figure D1 .
Figure D1.Illustration of the routine used to determine the height of largest crown extension at equilibrium (hlceeq) of a target tree in three contrasted situations of competition.A first step consists in determining the intersection between the potential crown of the target tree and the competitor.Then, the hlceeq is fixed between the maximum hlce (corresponding to the intersection between potential crowns) and the minimum hlce (which the height to crown base) based on the relative height of the competitor.5

Figure E1 .
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.5 hlce% fraction of the total height corresponding to the height of largest crown extension m determined based on tree inventory data of the study site hcb% fraction of the total height corresponding to the crown base height m determined based on tree inventory data of the study site Dd parameters of the crown to stem diameter function (Eq.10) m m -116.20/0.0280/0.00/0.0010.49/0.00/1379/-2881determined based on tree inventory data of the study site sh coefficient used to shift the mean crown to stem diameter ratio relationship to its maximum dimensionless 1

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 Estimation of hlce at equilibrium for a competitor of higher size Estimation of hlce at equilibrium for a competitor of lower size