HETEROFOR 1.0: a spatially explicit model for exploring the response of structurally complex forests to uncertain future conditions – Part 2: Phenology and water cycle

Climate change affects forest growth in numerous and sometimes opposite ways, and the resulting trend is often difficult to predict for a given site. Integrating and structuring the knowledge gained from the monitoring and experimental studies into process-based models is an interesting approach to predict the response of forest ecosystems to climate change. While the first generation of models operates at stand level, one now needs spatially explicit individual-based approaches in order to account for individual variability, local environment modification and tree adaptive behaviour in mixed and uneven-aged forests that are supposed to be more resilient under stressful conditions. The local environment of a tree is strongly influenced by the neighbouring trees, which modify the resource level through positive and negative interactions with the target tree. Among other things, drought stress and vegetation period length vary with tree size and crown position within the canopy.

In this paper, we describe the phenology and water balance modules integrated in the tree growth model HETEROFOR (HETEROgenous FORest) and evaluate them on six heterogeneous sessile oak and European beech stands with different levels of mixing and development stages and installed on various soil types. More precisely, we assess the ability of the model to reproduce key phenological processes (budburst, leaf development, yellowing and fall) as well as water fluxes.

Two two-phase models differing regarding their response function to temperature during the chilling period (optimum and sigmoid functions) and a simplified one-phase model are used to predict budburst date. The two-phase model with the optimum function is the least biased (overestimation of 2.46 d), while the one-phase model best accounts for the interannual variability (Pearson's r=0.68). For the leaf development, yellowing and fall, predictions and observations are in accordance. Regarding the water balance module, the predicted throughfall is also in close agreement with the measurements (Pearson's r=0.856; bias $=-\mathrm{1.3}$ %), and the soil water dynamics across the year are well reproduced for all the study sites (Pearson's r was between 0.893 and 0.950, and bias was between −1.81 and −9.33 %). The model also reproduced well the individual transpiration for sessile oak and European beech, with similar performances at the tree and stand scale (Pearson's r of 0.84–0.85 for sessile oak and 0.88–0.89 for European beech). The good results of the model assessment will allow us to use it reliably in projection studies to evaluate the impact of climate change on tree growth in structurally complex stands and test various management strategies to improve forest resilience.

1 Introduction

Climate projections for the future indicate a substantial increase in air temperature all over Europe (between 1.0 and 5.5 ^∘C, depending on the greenhouse gas emission scenario) and changes in precipitation regime according to the region (Jacob et al., 2014; Kovats et al., 2014). Climate extremes (e.g. heat waves and droughts) are also predicted to increase in intensity and frequency (Dai, 2013; Jacob et al., 2018). These changing climate conditions affect forest growth and mortality (Allen et al., 2015; Teskey et al., 2015; Charru et al., 2017; Kornhuber et al., 2019) and have an impact on the provision of ecosystem services (Hassan et al., 2005; Shvidenko et al., 2005; Rasche et al., 2013). Among these services, forests play an important role in regulating the climate system by sequestering carbon in biomass and soil (Myhre et al., 2013; Le Quéré et al., 2018) and by determining water and energy exchanges with the atmosphere through their evapotranspiration and land surface properties (e.g. albedo, roughness) (Bonan, 2008; Stocker et al., 2013).

Since climate change affects some tree growth processes positively and others negatively and given the interactions among factors as well as the feedback and acclimation mechanisms, it is not easy to predict the resulting effect on tree growth at a given site (Lindner et al., 2014; Herr et al., 2016). Knowledge about climate change has been acquired based on long-term monitoring studies that are limited to the observed changes (Bussotti and Pollastrini, 2017; Etzold et al., 2019) and on experiments of environment manipulation generally analysing one or two factors at a time for a limited period (Ainsworth and Long, 2005; Norby et al., 2010; Wolkovich et al., 2012; Meir et al., 2015). In order to apprehend the complex functioning of forest ecosystems, the use of process-based modelling is a complementary approach that allows the integration and structuring of the existing knowledge and extrapolations to be made for unprecedented conditions like those projected for the coming decades.

Process-based models were originally built to predict forest growth response to environmental changes at stand level without accounting for management operations and canopy heterogeneity. Such models were therefore suitable for pure even-aged stands but hardly manage to simulate mixed and structurally complex stands (Dufrêne et al., 2005; Pretzsch et al., 2007). Yet, nowadays, a promising way to adapt forests to climate change is to progressively turn them into uneven-aged and mixed stands using continuous cover forestry and natural-disturbance-based management to improve their stress resistance and resilience (DeRose and Long, 2014; Messier et al., 2015; Anderegg et al., 2018). To account for the spatial heterogeneity, some process-based models were designed or adapted to simulate various tree cohorts (Collalti et al., 2016). However, this approach only considers the vertical dimension of spatial heterogeneity while implementing innovative forestry practices in structurally complex stands requires the horizontal dimension to be accounted for through a spatially explicit approach at tree level (Pacala and Deutschman, 1995; Pretzsch et al., 2007; Berger et al., 2008; Bravo et al., 2019).

To reproduce the complexity of forest ecosystem functioning in mixed and structured forests, models must take individual variability, local environment and tree adaptive behaviour into account (Berger et al., 2008). Tree size and species influence physiological and morphological properties that in turn affect the main growth processes (Binkley et al., 2013). Considering average individuals is therefore a rough approximation and does not allow for all the variability within a heterogeneous forest to be accounted (Berger et al., 2008). Even in cohort-based approaches, tree grouping can only be done for a limited number of criteria that are not necessarily representative of the whole tree diversity. The local environment of a tree is strongly influenced by the neighbouring trees, which modify the resource level through positive and negative interactions with the target tree (Grossiord et al., 2014). As trees compete for limited resources, neighbouring trees can decrease light, water and nutrient availability. Tree species can however develop strategies to avoid competition by using different temporal and spatial niches (complementarity; Grossiord, 2018). Positive interactions may also occur when the neighbouring trees improve the growing conditions of the target trees (facilitation; Pretzsch et al., 2013). Finally, trees adapt their morphology and physiological behaviour to the local environmental conditions by optimizing carbon allocation in order to maximize the acquisition of the limiting resource (Petritan et al., 2009; Yuang et al., 2019).

As this study focus on phenology and water cycling, we briefly review how these processes are influenced by tree characteristics and local environment. Phenology timing varies among tree species, which favours early-leafing species but can also expose them to late frosts (Lopez et al., 2008; Liu et al., 2018). Many studies report that leaf development starts earlier and leaf senescence occurs later in the understory compared to the overstorey (Gill et al., 1998; Seiwa, 1999a; Augspurger and Bartlett, 2003; Schieber, 2006; Vitasse, 2013; Gressler et al., 2015), which allows the understory trees to benefit from a longer growing period and consequently, to increase their productivity (Jolly et al., 2004). The presence of warmer temperatures in the understory is one of the hypotheses advanced to explain this difference in budburst between under- and overstorey (Augspurger and Bartlett, 2003; Schieber, 2006). Using a construction crane, Vitasse (2013) tested this hypothesis by transplanting seedlings of five tree species at 30 and 35 m height in the canopy. He observed that the budburst of the seedling growing at these heights was much earlier than that of the dominant trees. He concluded that the main factor to explain this difference in budburst is driven by ontogeny (tree age and height) as stated by Seiwa (1999b) and that the vertical profile in temperature within the canopy only plays a secondary role. To capture the differences in budburst between understorey and dominant trees, ontogeny must be taken into account in priority.

Drought stress occurs when trees can no longer adjust their water use to soil water availability, which reduces growth and can even lead to mortality in the short or medium term due to hydraulic failure or progressive carbon starvation (McDowell and Allen, 2015; Meir et al., 2015; Greenwood et al., 2017). The stomatal control of water use varies among tree species and depends on tree size (Martínez-Vilalta and Lloret, 2016). In general, stomatal conductance decreases with tree height, which can be related to the fact that taller trees experience higher hydraulic resistance, higher soil-to-leaf water potential differences and are more vulnerable to cavitation (Grote et al., 2016). For the same climate conditions above the forest canopy, water demand varies with the degree of crown shading (local microclimate), which depends on the crown position within the canopy (Bennett et al., 2015). All in all, dominant trees are more susceptible to drought stress and mortality since they are more exposed to stressful conditions (excessive radiation, high vapour pressure deficit and elevated temperature) and present a higher risk of cavitation (Grote et al., 2016; Rötzer et al., 2017). In addition, as dominant trees have higher evapotranspiration rates, the soil water reserves in their surroundings are more rapidly depleted, which is however partly compensated by deeper rooting and horizontal water redistribution. These dominant trees reduce water availability for suppressed trees but, at the same time, decrease their evaporation demand. Complementarity in water use can occur when trees of different size and species take up water from different soil layers (Schwendenmann et al., 2015). This can also result in facilitation through hydraulic lift (Zapater et al., 2011). Mixed and structured stands promote facilitation and complementarity in water use but can also lead to faster exploitation of soil water reserves (Schäfer et al., 2018).

Modelling the complex functioning of heterogeneous forests is rather challenging. A more detailed representation of tree interactions comes at the price of a higher complexity, eventually lower robustness and longer computing times. One needs however spatially explicit individual-based models for gaining a mechanistic and comprehensive understanding of tree interactions and for comparing various spatial representations of stand structure in order to select the best one for the considered function (Berger et al., 2008; Bravo et al., 2019). Among other things, such models allow tree spatial configuration to be taken into account and stands composed of the same trees but with a contrasted spatial aggregation to be distinguished between (e.g. intimate vs. patch-wise mixture). However, very few models of this type currently exist. For all of these reasons, we decided therefore to develop a spatially explicit individual-based model called HETEROFOR for HETEROgeneous FOrest.

The processes regulating the carbon fluxes and the dimensional growth constitute the core of the HETEROFOR model and are described in Jonard et al. (2020a). Here, we focus on the description of two modules essential for predicting the impact of climate change on tree growth: phenology and water balance (Park et al., 2016; Choat et al., 2018). Phenology is described at the species level, with the possibility to make it dependent on tree size. Water balance can be achieved at the tree level or at the stand level by aggregation of individual tree properties. We used data from long-term forest monitoring to evaluate the capacity of the model to reproduce key phenological phases (budburst, leaf development, yellowing and fall) and the soil water content dynamics, as well as to estimate individual transpiration, stand throughfall and deep drainage. Evaluating each module separately is necessary to ensure the consistency of the whole model (Soares et al., 1995).

2 Material and methods

2.1 Model description

2.1.1 Overall model

HETEROFOR is a model hosted in Capsis (Computer-Aided Projections of Strategies In Silviculture), a software platform for forest growth simulations (Dufour-Kowalski et al., 2012) that provides the execution system and procedures to run simulations and display the outputs. Still, apart from these data structures and operative methods, all initialization and evolution procedures are specific to HETEROFOR. The initialization phase of the model consists in loading different files (tree species parameters, tree and stand characteristics, chemical and physical soil properties, meteorological data, and fruit production data) in order to create trees and soil horizons. Then, tree growth is calculated yearly according to the HETEROFOR methods presented in Jonard et al. (2020a). So far, HETEROFOR is adapted and calibrated only for deciduous species, but the adaptation to evergreen species is under progress.

Once the initialization is completed, the first routine called is the calculation of phenological periods from meteorological data, which is described is Sect. 2.1.2. This function provides key phenological dates and the daily foliage state (foliage development stage and green vs. discoloured leaf proportion) during the year. These phenological outputs are notably used for the radiation budget carried out using the SamsaraLight library coupled to HETEROFOR (Courbaud et al., 2003). According to a ray tracing approach and based on the solar radiation from the meteorological file, this library differentiates the direct and the diffuse components of the global radiation and determines, for both, the part of energy absorbed by the crown and the trunk of each tree and the part transmitted to the forest floor. The intercepted radiation is required to estimate evapotranspiration and tree photosynthesis. All aboveground and belowground water fluxes are calculated according to the processes described in Sect. 2.1.3, which allows the performance of an hourly water balance for each soil horizon at the tree or stand scale.

For each tree, gross primary production (gpp) is estimated either annually with a radiation use efficiency approach or daily using the photosynthesis method implemented in the model CASTANEA of Capsis (Dufrêne et al., 2005). In the latter case, the daily gpp is cumulated over the year. At the end of the year, a part of the annual gpp is used for growth and maintenance respiration, the remaining part constituting the net primary production (npp). Maintenance respiration can be estimated as a fraction of the gpp or calculated for each tree compartment by a method accounting for the living biomass, its nitrogen concentration and a Q₁₀ function that describes the temperature dependence. Growth respiration corresponds to a fraction of the carbon used to build the new tissues. The npp is then distributed to the different tree compartments (branches, trunk, roots and leaves), giving priority to the functional organs, namely leaves, fine roots and fruits. The carbon sharing between leaves and fine roots depends on the tree nutritional status, trees with a poorer nutrient status allocating relatively more carbon to fine roots. After carbon allocation to leaves, fine roots and fruits, the residual npp is distributed to structural tree parts (stem, branches and coarse roots) based on biomass allometry relationships. All of these processes involving carbon fluxes are described in detail in Jonard et al. (2020a). The HETEROFOR model also contains a tree nutrition and nutrient cycling module that will be described later.

2.1.2 Phenological module

The phenological module aims at predicting the temporal variation in the foliage status during the vegetation period. From budburst, leaf biomass progressively increases until a maximum value, then remains constant and finally decreases during leaf fall. This temporal evolution is characterized by the proportion of leaf biomass relative to its maximum value at full leaf development. In addition, two types of leaves are distinguished: green and discoloured leaves. The green leaf proportion is the ratio between the green leaf and the maximum leaf biomass. These two foliage properties are key variables used to simulate energy, water and carbon fluxes within the forest ecosystem. Photosynthesis and tree transpiration are dependent on the proportion of green leaves, since they are not active anymore on discoloured leaves. When leaves start yellowing, they still intercept rainfall, while their photosynthetic activity and transpiration are progressively reduced.

The following phenological phases are distinguished, in chronological order:

• chilling period or endodormancy: accumulation of coldness that breaks the bud dormancy; it is initiated at the chilling starting date (t₀) and ends at the forcing starting date (t₁);

• forcing period or ecodormancy: accumulation of heat that initiates the leaf development in the bud and leads to the budburst (budburst date is t_2a);

• leaf development: progressive growth of the leaves from budburst to the complete leaf development (leaf development date is t_2b);

• ageing: accumulation of coldness that is initiated at the ageing starting date (t₃) and ends at the yellowing starting date (t_4a);

• yellowing: loss of photosynthetic activity linked to the decrease in day length; this phase ends at the yellowing ending date (t_4b); and

• falling: the fall of the dead leaves starts (t_5a) when less than 60 % of the leaves are still green and continues until the leaf fall ending date (t_5b).

Since the phenological timing can vary considerably between species, the phenology dates are calculated for each tree species separately. Intra-specific differences are also likely to occur according to the size or social status (Cole and Sheldon, 2017) and can be optionally accounted for as described later.

The phenological module is optional in HETEROFOR. Activating the phenology requires an hourly meteorological file. If not activated, the model uses the budburst and leaf fall dates provided by the user, which are identical for all years and tree species.

The principle behind the whole phenology module is similar for each phase. A state variable is increasing progressively growing at a rate depending on meteorological conditions (air temperature). When the phase state reaches a certain threshold, the start of a new phase is triggered, except for the leaf yellowing and fall, which are partly simultaneous.

A two-phase model considering chilling and forcing is implemented to calculate the average budburst date (t_2a). The model starts to operate when the day of the year corresponds to the chilling starting date (t₀). At this moment, the daily chilling rate (R_c) is calculated according to a response function to temperature. Two different response functions can be chosen. The first, called optimum, considers minimum, maximum and optimal temperatures regarding the accumulation of coldness according to the following equation based on Hänninen (1990) as follows:

$$\begin{array}{}\text{(1)}& R_{\mathrm{c}}=\left\{\begin{array}{ll}\mathrm{0},& T\le T_{\min }\\ \frac{T-T_{\min }}{T_{\mathrm{opt}}-T_{\min }},& T_{\min }<T\le T_{\mathrm{opt}}\\ \frac{T-T_{\max }}{T_{\mathrm{opt}}-T_{\max }},& T_{\mathrm{opt}}<T\le T_{\max }\\ \mathrm{0},& T\ge T_{\max },\end{array}\right.\end{array}$$

with T_min, T_max and T_opt, the minimum, maximum and optimal temperatures (in degrees Celsius), respectively, and T, the daily average temperature (in degrees Celsius).

The second response function to temperature uses a sigmoid function (Chuine, 2000) as follows:

$$\begin{array}{}\text{(2)}& R_{\mathrm{c}}=\left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{1}+e^{C_a(T-C_c{)}^{\mathrm{2}}+C_b(T-C_c)}},& -\mathrm{5}\le T\le \mathrm{10}\\ \mathrm{0},& T>\mathrm{10}\mathrm{or}T<-\mathrm{5},\end{array}\right.\end{array}$$

with C_a, C_b and C_c (in degrees Celsius), chilling parameters.

This rate is summed each day until reaching the chilling threshold (C^*) that triggers the forcing process and sets the forcing starting date (t₁) to the current day. Regarding the forcing period, the forcing rate (R_f) is calculated using another sigmoid equation (Chuine, 2000) as follows:

$$\begin{array}{}\text{(3)}& R_{\mathrm{f}}=\left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{1}+e^{F_b(T-F_c)}},& T>T_{\mathrm{b}\mathrm{\_}\mathrm{for}}\\ \mathrm{0},& T\le T_{\mathrm{b}\mathrm{\_}\mathrm{for}},\end{array}\right.\end{array}$$

with F_b and F_c (in degrees Celsius), forcing parameters, and T_{b_for}, the base temperature for forcing.

The budburst is activated when the sum of the daily forcing rates reaches the forcing threshold (F^*).

A simplified one-phase version is implemented as well, which only considers forcing similarly to the two-phase model (Eq. 3). In this case, the forcing starting date (t₁) must be provided.

As the module was calibrated based on observations carried out on trees representative of the stand, the predicted budburst starting date is expected to be that of an average tree. Since, at this date, the leaf expansion of some trees has already started in real conditions, the model shifts the budburst date to correspond to that of the earliest trees. This budburst shift, t2a_shift, is equal to half the period between the budburst of the first and the last tree and must be provided by the user for the various tree species. By doing so, leaf development starts early for all trees which follow a same average evolution when belonging to a same-tree species.

Once the budburst starting date (t_2a) is calculated, the equations for the subsequent phenological variables are the same. The leaf development rate (R_ld) is cumulated daily until the leaf development threshold (LD^*) is reached. It is computed according to

$$\begin{array}{}\text{(4)}& R_{\mathrm{ld}}=\left\{\begin{array}{ll}T,& T>\mathrm{0}\\ \mathrm{0},& T\le \mathrm{0},\end{array}\right.\end{array}$$

where T is the daily average temperature of the current day (in degrees Celsius).

The leaf proportion (leafProp; in grams per gram) is calculated daily for each tree species (sp) according to

$$\begin{array}{}\text{(5)}& {\mathrm{leafProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}}={\displaystyle \frac{\sum _{t_{\mathrm{2}\mathrm{a}}}^t{R}_{\mathrm{ld}}}{{\mathrm{LD}}^{*}}},\end{array}$$

with t, the current day.

As many studies have shown that budburst in the understory occurs earlier than in the overstorey and ascribed this primarily to ontogeny (Gill et al., 1998; Seiwa 1999a, b; Augspurger and Bartlett, 2003; Schieber, 2006; Vitasse, 2013), we implemented an option to make the phenology size-dependent (phenology at tree level). With this option, the leaf development is first triggered in the smallest trees of each tree species and then progressively in the tallest ones according to their height. At the stand level, the option phenology at tree level provides exactly the same leaf development as the default option, but the differences appear at the tree scale. The default option assumes that all trees of a same species initiate budburst at the same time and display the same progressive leaf development, while the alternative one supposes that trees break down one after the other depending on their size.

With the option phenology at tree level, the leaf proportion of each tree (leafProp_{tree_t}) is updated daily (t) between the budburst starting date (t_2a) and the budburst ending date (t_2b) based on the leaf proportion calculated at the stand scale for the corresponding tree species (leafProp_{sp_t}) by

$$\begin{array}{}\text{(6)}& {\mathrm{leafProp}}_{\mathrm{tree}\mathrm{\_}\mathrm{t}}=\left\{\begin{array}{ll}\mathrm{1},& \frac{{\sum }_{\mathrm{1}}^{\mathrm{tree}}{a}_{\mathrm{leaf}}}{A_{\mathrm{leaf}}}\le {\mathrm{leafProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}}\\ \mathrm{0},& \frac{{\sum }_{\mathrm{1}}^{\mathrm{tree}}{a}_{\mathrm{leaf}}}{A_{\mathrm{leaf}}}>{\mathrm{leafProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}},\end{array}\right.\end{array}$$

with “tree”, the tree of interest (note that the trees are sorted in ascending order based on their height), a_leaf, the tree leaf area (in square metres), and A_leaf, the total stand leaf area (in square metres).

A fixed date, defined according to Dufrêne et al. (2005), is considered for the start of the ageing process (t₃). This process does not alter leaf quality but is a prerequisite for leaf yellowing (t_4a) that is initiated when the cumulated daily ageing rate (R_age) equals the ageing threshold (A^*).

$$\begin{array}{}\text{(7)}& R_{\mathrm{age}}=\left\{\begin{array}{ll}{T}_{\mathrm{b}\mathrm{\_}\mathrm{age}}-T,& T<T_{\mathrm{b}\mathrm{\_}\mathrm{age}}\\ \mathrm{0},& T\ge T_{\mathrm{b}\mathrm{\_}\mathrm{age}},\end{array}\right.\end{array}$$

with T_{b_age}, the base temperature for ageing (in degrees Celsius).

The leaf yellowing calculation gives the green leaf proportion, greenProp (in grams per gram), which provides the fraction of remaining green leaves compared to the maximum amount of green leaves for each tree species. It is set to 1 before the start of yellowing and then decreases with day length according to the following equation:

$$\begin{array}{}\text{(8)}& {\mathrm{greenProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}}={\mathrm{greenProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}-\mathrm{1}}\times {\left({\displaystyle \frac{{\mathrm{DL}}_t-{\mathrm{DL}}_{\min }}{{\mathrm{DL}}_{t_{\mathrm{4}\mathrm{a}}}-{\mathrm{DL}}_{\min }}}\right)}^y,\end{array}$$

with DL_t and DL${}_{t_{\mathrm{4}\mathrm{a}}}$, the day lengths (in hours) for the current day and t_4a, respectively, DL_min,the minimum day length (in hours) value over the year, and y, a leaf yellowing parameter.

The day length (in hours) is calculated according to Teh (2006):

$$\begin{array}{}\text{(9)}& \mathrm{DL}={\displaystyle \frac{\mathrm{24}}{\mathit{\pi }}}\times \mathrm{acos}\left(-{\displaystyle \frac{\sin \left(\mathit{\delta }\right)\times \sin \left(\mathit{\lambda }\right)}{\cos \left(\mathit{\delta }\right)\times \cos \left(\mathit{\lambda }\right)}}\right),\end{array}$$

where λ is the site latitude (in radians) and δ, the solar declination (in radians) determined as $\mathit{\delta }=-\frac{\mathrm{23.45}\times \mathit{\pi }}{\mathrm{180}}\times \cos \left(\mathrm{2}\mathit{\pi }\frac{\mathrm{DOY}+\mathrm{10}}{\mathrm{365}}\right)$ and DOY is the day of the year (e.g. 1 January is 1, 2 January is 2 and 1 February is 32).

The yellowing phase ends when the green leaf proportion drops below a threshold, called yellowing threshold, Y^*, indicated by the model user in the species file. The leaf fall (t₅) is set to start rapidly after yellowing initiation, namely when greenProp reaches 0.60, considering that leaves that are no longer photosynthetically active can quickly fall.

The falling rate (R_fall) is calculated daily and is used to update leafProp for each tree species. It depends on the wind and frost episodes. While the frost weakens the leaf petiole, the wind can break it and take away the leaf. For this reason, leafProp is determined for each day t as follows:

$$\begin{array}{}\text{(10)}& {\mathrm{leafProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}}={\mathrm{leafProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}-\mathrm{1}}-f_{\mathrm{ampl}}\times \mathrm{WS}\times R_{\mathrm{fall}},\end{array}$$

with f_ampl, a frost amplifier coefficient fixed to 1 before the occurrence of 5 consecutive hours with air temperature below 0 ^∘C and is then set to 2 and 3 for oak and beech, respectively; WS is the daily average wind speed (in metres per second), and R_fall is the falling rate (in seconds per metre per day) calibrated as described in Sect. 2.2.

According to Eq. (10), leafProp_{sp_t} progressively decreases from 1 to 0, but it cannot take a value below greenProp_{sp_t}, accounting for the fact that green leaves are not expected to fall. Finally, when all leaves have fallen, the trees enter in the leafless period until the budburst of the following year.

Similarly to leaf development but with a reverse order, the option phenology at tree level first triggers the leaf yellowing and fall in the taller trees and then in the smaller ones in order to reproduce the observations reported by Gressler et al. (2015). This options updates daily the green leaf and leaf proportions of each tree (greenProp_{tree_t}, leafProp_{tree_t}) between the yellowing starting date (t_4a) and the falling ending date (t_5b) based on the green leaf and leaf proportions calculated at the stand scale for the corresponding tree species (greenProp_{sp_t}, leafProp_{sp_t}) as follows:

$$\begin{array}{}\text{(11)}& {\mathrm{greenProp}}_{\mathrm{tree}\mathrm{\_}\mathrm{t}}=\left\{\begin{array}{ll}\mathrm{1},& \frac{{\sum }_n^{\mathrm{tree}}{a}_{\mathrm{leaf}}}{A_{\mathrm{leaf}}}\le {\mathrm{greenProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}}\\ \mathrm{0},& \frac{{\sum }_n^{\mathrm{tree}}{a}_{\mathrm{leaf}}}{A_{\mathrm{leaf}}}>{\mathrm{greenProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}}\end{array}\right.\end{array}$$

$$\begin{array}{}\text{(12)}& {\mathrm{leafProp}}_{\mathrm{tree}\mathrm{\_}\mathrm{t}}=\left\{\begin{array}{ll}\mathrm{1},& \frac{{\sum }_n^{\mathrm{tree}}{a}_{\mathrm{leaf}}}{A_{\mathrm{leaf}}}\le {\mathrm{leafProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}}\\ \mathrm{0},& \frac{{\sum }_n^{\mathrm{tree}}{a}_{\mathrm{leaf}}}{A_{\mathrm{leaf}}}>{\mathrm{leafProp}}_{\mathrm{sp}\mathrm{\_}\mathrm{t}},\end{array}\right.\end{array}$$

with tree, the tree of interest (note that the trees are sorted by descending order based on their height), a_leaf, the tree leaf area (in square metres), and A_leaf, the total stand leaf area (in square metres).

The option phenology at tree level provides the opportunity to compare two contrasted hypotheses regarding individual tree phenology and to evaluate to what extent it has an impact on tree growth.

2.1.3 Water balance module

The water balance module operates at an hourly time step and simulates the partitioning of incident rainfall into the main forest water fluxes and pools, namely interception (i.e. water storage on foliage and bark and evaporation), throughfall, stemflow, water movements between soil horizons and deep drainage, transpiration and soil water uptake in the different soil horizons, and soil evaporation (Fig. 1). Surface runoff and groundwater level rise are not yet considered in the current HETEROFOR version. Instead, when saturation is reached in a soil layer, the water surplus is transferred to the horizon below or is lost when it occurs in the last horizon.

Figure 1 Schematic representation of the water fluxes and pools in the water balance module. Rainfall is divided into throughfall directly reaching the forest floor and a pre-stemflow component intercepted by the foliage and the bark. Once the foliage and bark are saturated, the water surplus increases the throughfall flux and flows along the branches and the trunk to generate stemflow. The throughfall and stemflow fluxes enter in the upper part of the soil and then, move from one horizon to the other according to the Darcy's law. For a soil horizon, hr, the water input fluxes can be the drainage from the upper horizon (D_hr−1) and the capillary rise from the lower horizon (CR_hr+1) that depend on the water potential gradient between the concerned horizons and on their hydraulic conductivity. The output fluxes are the drainage (D_hr) and the capillary rise (CR_hr), the root water uptake (UP_root(hr)), and the surplus (S_hr) that appears when the horizon water content exceeds the saturated water content. One part of this last flux can directly leave the system as deep drainage (DD) when preferential flow is considered, in addition to the water drainage of the last horizon. In parallel, water evaporates from foliage, bark and soil and is taken up by roots to enable tree transpiration. The evapotranspiration fluxes are all calculated with the Penman–Monteith equation.

In a first step, the parameters considered as constant during the leaved and leafless periods are estimated. Then, the various water fluxes are calculated at an hourly time step. The default option for the water balance module calculates the water fluxes at the stand level by summing properties estimated at the tree level (maximum foliage and bark storage capacities and throughfall and stemflow proportions). For this option, tree transpiration is calculated at the tree level and summed at the stand scale. Stand transpiration is then used to estimate root water uptake in the different soil horizons, assuming that all trees are taking up water in the same reservoirs in which soil water is redistributed homogeneously between two hourly time steps. This hypothesis can be justified by soil anisotropy, which induces a higher horizontal than vertical soil conductance. This is justified since water movements through the same horizon depend only on its own hydrological properties, while the presence of one horizon with a low conductance can slow down vertical water movement in the upper horizons (Todd and Mays, 2005). Moreover, as sediments are preferentially deposited on their longest side, the vertical conductance is decreased with regard to the horizontal one (Cristiano et al., 2016) so that the ratio of the horizontal vs. vertical conductance ranges between 2 and 10 in alluvial soils and amounts to 100 in clay soils (Todd and Mays, 2005).

The user can select an alternative option activate fine spatial resolution to perform water balance on an individual scale. In this case, all the water fluxes (throughfall, stemflow, foliage, bark and soil evaporation, transpiration, water uptake, soil water movements and drainage) are calculated at the individual level. For this option, the model distributes the total soil volume in individual soil volumes (called pedons).

The pedon area (a_pedon) is determined proportionally to the leaf area of the associated tree (but is limited to 2 times its crown projection) as follows:

$$\begin{array}{}\text{(13)}& a_{\mathrm{pedon}}={\displaystyle \frac{a_{\mathrm{leaf}}}{A_{\mathrm{leaf}}}}{A}_{\mathrm{stand}},\end{array}$$

with a_leaf, the tree leaf area (in square metres), A_leaf, the total stand leaf area (in square metres), and A_stand, the total stand area (in square metres).

In sparse stands, all the stand area is not allocated to the trees, and the remaining area is considered as a pedon without any associated tree. With the fine spatial resolution, the model performs a water balance for each tree pedon and also for the remaining pedons (without tree). Contrary to the default option, the alternative option supposes no water redistribution among pedons. This hypothesis could become more appropriate than the perfect-redistribution hypothesis when soil dries (Friedman and Jones, 2001), at least beyond the air entry value (Assouline and Or, 2006). The two options allow the user to test two contrasting hypotheses regarding soil water redistribution in the horizontal dimension. In the following description, variables calculated at the stand scale are represented with capital letters, while lowercase letters are used for variables at the tree level. In some cases, when the equation is the same at the tree and the stand level, the variables are represented only with capital letters to avoid unnecessary duplications.

Foliage and bark storage capacity

The maximum foliage storage capacity of a tree (c_{foliage_max}; in litres) is calculated by multiplying the foliage storage capacity of the corresponding tree species by the tree leaf area as follows:

$$\begin{array}{}\text{(14)}& c_{\mathrm{foliage}\mathrm{\_}\max }=a_{\mathrm{leaf}}\cdot c_{\mathrm{foliage}\mathrm{\_}\mathrm{sp}},\end{array}$$

with c_{foliage_sp}, the foliage storage capacity for the species sp (in millimetres or litres per square metre of leaves).

To obtain it at the stand level (C_{foliage_max}; in litres), the model sums the maximum foliage storage capacity of all the trees.

Bark storage capacity depends on season (i.e. leaved and leafless periods) and on tree species. It is derived from a linear model proposed by André et al. (2008a) predicting the individual stemflow (sf; in litres) produced during a rain event as a function of tree girth (C130; in centimetres) and rainfall amount (R; in millimetres).

$$\begin{array}{}\text{(15)}& \mathrm{sf}=a+b\cdot \mathrm{C}\mathrm{130}+c\cdot R+d\cdot \mathrm{C}\mathrm{130}\cdot R+\mathit{\tau }+\mathit{\delta }+\mathit{\epsilon },\end{array}$$

where a (in litres), b (in litres per centimetre), c (in square metres) and d (in square metres per centimetres) are fixed effect parameters varying with tree species and season; τ and δ are random factors characterizing the tree and rain event variability, respectively, and ε account for the residuals.

As it is multiplied by the rainfall amount in Eq. (15), the term “$c+d\cdot \mathrm{C}\mathrm{130}$” may be interpreted as an estimate of the stemflow rate (sf_rate; in litres per millimetre). In parallel, André et al. (2008a) determined the rainfall threshold for stemflow appearance (R_min; in millimetres), defined as the amount of rainfall required to produce stemflow at the base of the trunk. This threshold was found to be independent of tree size, while it depends on both season and tree species. Multiplying the sf_rate estimations by R_min values for the corresponding species and season provides estimates of the tree bark storage capacity (c_bark; in litres), namely the amount of water accumulated on branch and trunk bark before stemflow occurs at tree base, calculated as follows:

$$\begin{array}{}\text{(16)}& c_{\mathrm{bark}}=(c+d\cdot \mathrm{C}\mathrm{130})\cdot R_{\min }.\end{array}$$

The individual c_bark estimates are then summed over all trees of a same species for each season to determine leafless (ll) and leaved (ld) stand bark storage capacities (C_{bark_sp_ll} and C_{bark_sp_ld}; in litres). As shown by André et al. (2008a), the seasonal variation in the bark storage capacity is not significant since the corresponding changes in the three parameters (c, d and R_min) offset each other. We maintained however the distinction between seasons since the parameters of Eq. (15) were also used to estimate throughfall and stemflow proportions (described hereafter), which are clearly season-dependent.

Throughfall and stemflow proportions

For a given tree, the proportion of stand rainfall reaching the ground at the base of the trunk as stemflow may be calculated by dividing the stemflow rate (see above) by the pedon or stand area (a_pedon or A_stand; in square metres) depending on the selected option (tree vs. stand scale water balance) as follows:

$$\begin{array}{}\text{(17)}& \mathit{\%}\mathrm{sf}={\displaystyle \frac{c+d\cdot \mathrm{C}\mathrm{130}}{a_{\mathrm{pedon}}}}\mathrm{or}{\displaystyle \frac{c+d\cdot \mathrm{C}\mathrm{130}}{A_{\mathrm{stand}}}}.\end{array}$$

For the water balance at the stand scale, the stemflow proportion per tree species is then calculated separately for the leafless and the leaved periods (% SF_{sp_ll}, % SF_{sp_ld}) by summing the corresponding tree stemflow proportions. The stemflow proportion is also calculated at the stand scale for each period (% SF_ll, % SF_ld). Finally, tree and stand level throughfall proportions are obtained directly from the stemflow proportions:

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }\mathrm{\%}{\mathrm{tf}}_{\mathrm{ll}}=\mathrm{1}-\mathrm{\%}{\mathrm{sf}}_{\mathrm{ll}}\mathrm{or}\mathrm{\%}{\mathrm{TF}}_{\mathrm{ll}}=\mathrm{1}-\mathrm{\%}{\mathrm{SF}}_{\mathrm{ll}}\text{(19)}& {\displaystyle }& {\displaystyle }\mathrm{\%}{\mathrm{tf}}_{\mathrm{ld}}=\mathrm{1}-\mathrm{\%}{\mathrm{sf}}_{\mathrm{ld}}\mathrm{or}\mathrm{\%}{\mathrm{TF}}_{\mathrm{ld}}=\mathrm{1}-\mathrm{\%}{\mathrm{SF}}_{\mathrm{ld}}.\end{array}$$

Absorbed radiation proportions

During the leaved period, the radiation absorbed by the trees is provided by the SamsaraLight library for either the whole period (simplified radiation balance, the default option) or every hour of key phenological dates (detailed radiation balance, an alternative option). It may be determined by either considering absorption by tree crowns as a function of leaf area density and ray path length through the crown by applying the Beer–Lambert law or specifying relative crown radiation absorption coefficients for each species. At the tree scale, the proportion of incident radiation absorbed per unit of leaf area during the vegetation period ($\mathrm{\%}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{canopy}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}$) is calculated as the ratio of the radiation absorbed by the crown over the whole vegetation period (aRAD_{tree_crown}; in megajoules) divided by the corresponding incident radiation (RAD; in megajoules per square metre) and the tree leaf area as follows:

$$\begin{array}{}\text{(20)}& \mathrm{\%}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}={\displaystyle \frac{{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{crown}}}{\mathrm{RAD}\cdot a_{\mathrm{leaf}}}}.\end{array}$$

At the stand scale, this proportion is obtained by summing the radiation absorbed by each crown and dividing it by the incident radiation and the leaf area of the whole stand as follows:

$$\begin{array}{}\text{(21)}& \mathrm{\%}{\mathrm{aRAD}}_{\mathrm{stand}\mathrm{\_}\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}={\displaystyle \frac{\sum _{\mathrm{tree}}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{crown}}}{\mathrm{RAD}\cdot A_{\mathrm{leaf}}}}.\end{array}$$

Similarly, the proportion of incident radiation absorbed per unit of bark area is obtained, at the tree and stand scales, respectively, by

$$\begin{array}{}\text{(22)}& \mathrm{\%}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{bark}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}={\displaystyle \frac{{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{trunk}}}{\mathrm{RAD}\cdot a_{\mathrm{bark}}}},\end{array}$$

with aRAD_{tree_trunk}, the radiation absorbed by the trunk of a given tree (in megajoules), a_bark, the tree bark area (in square metres), and

$$\begin{array}{}\text{(23)}& \mathit{\%}{\text{aRAD}}_{\mathrm{stand}\mathrm{\_}\mathrm{bark}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}={\displaystyle \frac{\sum _{\mathrm{tree}}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{trunk}}}{\mathrm{RAD}\cdot A_{\mathrm{bark}}}},\end{array}$$

with A_bark, the stand bark area (in square metres). At both scales (tree and stand), the proportion of incident radiation transmitted to the understorey is the transmitted radiation (transRAD; in megajoules per square metre), determined as the difference between the incident radiation and the radiation absorbed by the tree(s) divided by the incident radiation as follows:

$$\begin{array}{}\text{(24)}& \mathrm{\%}\mathrm{transRAD}={\displaystyle \frac{\mathrm{transRAD}}{\mathrm{RAD}}}.\end{array}$$

The radiation transmitted to the understory is then partitioned into the radiation intercepted by the ground vegetation and that reaching the soil by applying Beer–Lambert law, considering the ground vegetation leaf area index (described later in the section Ground vegetation transpiration and soil evaporation).

In the following sections, all of these proportions are used to estimate the hourly absorbed or transmitted radiations based on the hourly incident radiation.

For the leafless period, the proportions of incident radiation intercepted by the trunks and the branches and transmitted to the understory are obtained based on the Beer–Lambert law using the bark area index (i.e. bark surface divided by the stand or pedon area, BAI, in square metres per square metre) calculated from the bark biomass, density and thickness as follows:

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }\mathrm{\%}{\mathrm{aRAD}}_{\mathrm{bark}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}={\displaystyle \frac{\mathrm{1}-\exp (-k\cdot \mathrm{BAI})}{\mathrm{BAI}}}\text{(26)}& {\displaystyle }& {\displaystyle }\mathrm{\%}\mathrm{transRAD}={\displaystyle \frac{\exp (-k\cdot \mathrm{BAI})}{\mathrm{BAI}}}.\end{array}$$

Interception and evaporation of water stored on foliage and bark

Based on the preceding calculations, the water balance module starts updating the different water fluxes and pools for every hourly time step. First water evaporation from foliage and from bark is computed using the Penman–Monteith (P–M) equation (Monteith, 1965) at the tree or stand scale. The latent heat flux density is calculated as follows:

$$\begin{array}{}\text{(27)}& \mathit{\lambda }E={\displaystyle \frac{\mathrm{\Delta }R+\frac{\mathit{\rho }\cdot c_{\mathrm{p}}\cdot \mathrm{VPD}}{r_{\mathrm{a}}}}{\mathrm{\Delta }+\mathit{\gamma }\left(\frac{r_{\mathrm{a}}+r_{\mathrm{s}}}{r_{\mathrm{a}}}\right)}},\end{array}$$

where λE is the latent heat flux density (in watts per square metre), λ is the water latent heat of vaporization and is 2 454 000 J kg⁻¹ (Teh, 2006), γ is the psychometric constant and is 0.658 mbar K⁻¹ (Teh, 2006), Δ is the slope of the saturated vapour pressure curve (in millibars per kelvin),

$$\begin{array}{}\text{(28)}& \mathrm{\Delta }\approx {\displaystyle \frac{\mathrm{d}{e}_{\mathrm{s}}\left(T\right)}{\mathrm{d}T}}={\displaystyle \frac{\mathrm{25029.4}\cdot \exp \left[\frac{\mathrm{17.269}.T}{T+\mathrm{237.3}}\right]}{(T+\mathrm{237.3}{)}^{\mathrm{2}}}},\end{array}$$

ρ is the moist air density and is 1.209 kg m⁻³, c_p is the moist air specific heat capacity and is 1010 J kg⁻¹ K⁻¹, T is the air temperature (in degrees Celsius), R is the absorbed radiation per unit of leaf or bark area (in watts per square metre of leaves/bark), r_a is the aerodynamic resistance (in seconds per metre), the inverse of aerodynamic conductance, g_a,

$$\begin{array}{}\text{(29)}& r_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{g_{\mathrm{a}}}},\end{array}$$

r_s is the surface resistance (in seconds per metre), the inverse of surface conductance, g_s,

$$\begin{array}{}\text{(30)}& r_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{g_{\mathrm{s}}}},\end{array}$$

and VPD, the vapour pressure deficit (in millibars or hectopascals), calculated as follows, is based on the air temperature and the relative humidity,

$$\begin{array}{}\text{(31)}& \mathrm{VPD}=e_{\mathrm{s}}\left(T\right)-e_{\mathrm{r}},\end{array}$$

with e_s, the saturated vapour pressure (in millibars),

$$\begin{array}{}\text{(32)}& e_{\mathrm{s}}\left(T\right)=\mathrm{6.1078}\exp \left[{\displaystyle \frac{\mathrm{17.269}T}{T+\mathrm{237.3}}}\right],\end{array}$$

e_r, the air vapour pressure (in millibars), and

$$\begin{array}{}\text{(33)}& e_{\mathrm{r}}={\displaystyle \frac{\mathrm{RH}}{\mathrm{100}}}{e}_{\mathrm{s}}\left(T_{\mathrm{r}}\right),\end{array}$$

where RH is the relative humidity (in 10⁻² hectopascals per hectopascal).

The radiation absorbed hourly per unit of leaf area ($\mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}$; in watts per square metre) is obtained by multiplying the proportion of incident radiation absorbed per leaf area unit by the hourly incident radiation (h_RAD; in watts per square metre) as follows:

$$\begin{array}{}\text{(34)}& \mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}=\mathrm{\%}{\mathrm{aRAD}}_{\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}\cdot \mathrm{h}\mathrm{\_}\mathrm{RAD}.\end{array}$$

Similarly, the hourly absorbed radiation per unit of bark area ($\mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{bark}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}$; in watts per square metre) is obtained by multiplying the proportion of incident radiation absorbed by the bark by the hourly incident radiation:

$$\begin{array}{}\text{(35)}& \mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{bark}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}=\mathrm{\%}{\mathrm{aRAD}}_{\mathrm{bark}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}\cdot \mathrm{h}\mathrm{\_}\mathrm{RAD}.\end{array}$$

The aerodynamic resistance is defined as the inverse of the aerodynamic conductance, which represents the ease for a water vapour molecule to get away from its original location once it has been evaporated. Similarly, the surface resistance is the inverse of surface conductance that represents the ease for water molecules to migrate through the surface–air interface. The aerodynamic resistance depends mainly on wind speed and turbulence, while the surface resistance is a function of the water diffusivity through the surface.

According to Teh (2006) and depending on the scale considered (tree or stand), the mean canopy air resistance may be obtained by integrating the canopy air conductance (g_a; in metres per second) values estimated at 11 height levels between the mid-crown or mid-canopy height and the dominant height for the foliage and between half of the total or dominant height and the dominant height for the bark:

$$\begin{array}{}\text{(36)}& g_{\mathrm{a}}=\mathrm{0.006}\cdot \sqrt{{\displaystyle \frac{\mathrm{WS}}{l_{\mathrm{sp}}}}},\end{array}$$

with l_sp, the mean leaf width, and WS, the wind speed (in metres per second).

The mid-canopy height is determined as the mid-height between the dominant height of the stand (h_d; in metres), defined as the mean total height of the 100 biggest trees per hectare, and the canopy base height (h_cb; in metres), defined as the mean height to crown base of the 100 smallest trees per hectare. At the tree scale, the integration is done between the mid-crown height and the total height for the foliage and between half of the total height and the total height for the bark.

WS is estimated at the different heights (h; in metres) based on the dominant-height wind speed (${\mathrm{WS}}_{h_{\text{d}}}$; in metres per second) and on the wind speed attenuation coefficient (α):

$$\begin{array}{}\text{(37)}& \mathrm{WS}={\mathrm{WS}}_{h_{\text{d}}}\cdot e^{-\left[\mathit{\alpha }\cdot \left(\mathrm{1}-\frac{h}{h_{\text{d}}}\right)\right]},\end{array}$$

where WS${}_{h_{\text{d}}}$ is calculated according to Jetten (1996), based on the measured wind speed and its height of measurement.

$$\begin{array}{}\text{(38)}& \mathrm{WS}\left(h\right)=\mathrm{WS}\left(z_{\mathrm{m}}\right)\cdot {\displaystyle \frac{\ln \left[(z_{\mathrm{e}}-d_{\mathrm{m}})/z_{\mathrm{0}\mathrm{m}}\right]}{\ln \left[(z_{\mathrm{m}}-d_{\mathrm{m}})/z_{\mathrm{0}\mathrm{m}}\right]}}\cdot {\displaystyle \frac{\ln \left[(h-d_{\mathrm{f}})/z_{\mathrm{0}f}\right]}{\ln \left[(z_{\mathrm{e}}-d_{\mathrm{f}})/z_{\mathrm{0}f}\right]}},\end{array}$$

where h is the height at which wind speed is estimated (in this case the dominant height), z_e is the reference height (in metres) fixed to 50 m, z_m is the wind speed measurement height (2.5 m), d_m is the surface roughness height (in metres) of the meteorological station, fixed to 0.08 m, z_0m is the zero plane displacement (in metres) of the meteorological station, fixed to 0.015 m, d_f is the surface roughness height (in metres) of the forest and estimated as 0.75⋅h_d, and z_0f is the zero plane displacement (in metres) of the meteorological station, fixed to 0.1⋅h_d.

While no surface resistance is considered for the foliage evaporation (infinite conductance), the bark conductance (in metres per second) depends on the bark storage at the previous time step (prevS_{bark_sp}; in litres) and the bark storage capacity (C_{bark_sp}; in litres) according to

$$\begin{array}{}\text{(39)}& \begin{array}{rl}& g_{\mathrm{s}\mathrm{\_}\mathrm{bark}\mathrm{\_}\mathrm{sp}}=\\ & g_{\mathrm{s}\mathrm{\_}\mathrm{bark}\mathrm{\_}\min }+(g_{\mathrm{s}\mathrm{\_}\mathrm{bark}\mathrm{\_}\max }-g_{\mathrm{s}\mathrm{\_}\mathrm{bark}\mathrm{\_}\min })\cdot {\displaystyle \frac{{\mathrm{prevS}}_{\mathrm{bark}\mathrm{\_}\mathrm{sp}}}{C_{\mathrm{bark}\mathrm{\_}\mathrm{sp}}}}\end{array}.\end{array}$$

The latent heat flux density is then converted to hourly water evaporation (EV; in litres per hour per square metre of leaves):

$$\begin{array}{}\text{(40)}& {\mathrm{EV}}_{{\text{foliage or bark\_m}}^{\mathrm{2}}}={\displaystyle \frac{\frac{\mathit{\lambda }E}{\mathit{\lambda }}}{d_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}}\cdot \mathrm{1000}\cdot \mathrm{60}\cdot \mathrm{60},\end{array}$$

with E, the mass of water evaporated (in kilograms per square metre per second), and $d_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$, the water density (998 kg m⁻³).

Hourly tree or stand foliage evaporation (EV_{foliage_stand}; in litres per hour) is obtained by multiplying EV_foliage from Eq. (40) by the tree or stand leaf area as follows:

$$\begin{array}{}\text{(41)}& {\mathrm{EV}}_{\mathrm{foliage}}={\mathrm{EV}}_{\mathrm{foliage}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}\cdot \left(a_{\mathrm{leaf}}\mathrm{or}{A}_{\mathrm{leaf}}\right).\end{array}$$

Similarly, hourly evaporation from bark (EV_bark; in litres per hour) is determined separately for each tree or tree species by

$$\begin{array}{}\text{(42)}& {\mathrm{EV}}_{\mathrm{bark}}={\mathrm{EV}}_{\mathrm{bark}\mathrm{\_}\mathrm{sp}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}\cdot \left(a_{\mathrm{leaf}}\mathrm{or}{A}_{\mathrm{bark}\mathrm{\_}\mathrm{sp}}\right),\end{array}$$

where A_{bark_sp} is the bark area for the tree species sp (in square metres).

Evaporation from foliage and from bark cannot be larger than the corresponding amounts of water stored on these surfaces, namely S_foliage (in litres) and S_{bark_sp} (in litres) (see next section). Therefore, the following conditions are set:

$$\begin{array}{}\text{(43)}& {\displaystyle }& {\displaystyle }{\mathrm{EV}}_{\mathrm{foliage}}=min({\mathrm{EV}}_{\mathrm{foliage}},S_{\mathrm{foliage}})\text{(44)}& {\displaystyle }& {\displaystyle }{\mathrm{EV}}_{\mathrm{bark}}=min({\mathrm{EV}}_{\mathrm{bark}},S_{\mathrm{bark}\mathrm{\_}\mathrm{sp}}).\end{array}$$

Partitioning of rainfall into interception, throughfall and stemflow

Rainfall passing through the canopy can be intercepted by the foliage, the branches and the stems of the tree(s). These reservoirs saturate progressively, and the water then flows along the trunks to the tree base(s) to produce stemflow or drips from the canopy to the ground as throughfall. For some of the parameters (i.e. storage capacities, stemflow proportions) showing contrasting values depending on the season, the leaved and the leafless periods are distinguished to describe these processes. In addition, several intermediate state variables are considered, namely

• tree or stand rainfall

  $$\begin{array}{}\text{(45)}& \left(R_{\mathrm{tree}\mathrm{or}\mathrm{stand}},\mathrm{in}\mathrm{litres}\right)=R\cdot \left(A_{\mathrm{pedon}}\mathrm{or}{A}_{\mathrm{stand}}\right);\end{array}$$

• foliage storage (S_foliage; in litres) corresponding to the amount of water stored on the tree or stand foliage;

• previous stand foliage storage (prevS_foliage; in litres) being the tree or stand foliage storage at the previous time step;

• remaining foliage storage capacity (RemC_foliage; in litres), defined as

  $$\begin{array}{}\text{(46)}& {\mathrm{RemC}}_{\mathrm{foliage}}=C_{\mathrm{foliage}}-({\mathrm{prevS}}_{\mathrm{foliage}}-{\mathrm{EV}}_{\mathrm{foliage}}),\end{array}$$

  with non-intercepted rainfall (unintR; in litres).

For the leaved period, the foliage storage and the non-intercepted rainfall are updated at every time step, considering various cases as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{if}({\mathrm{RemC}}_{\mathrm{foliage}}>\mathrm{0})\mathit{\{}\\ {\displaystyle }& {\displaystyle }\mathrm{if}({\mathrm{RemC}}_{\mathrm{foliage}}>R_{\mathrm{tree}\mathrm{or}\mathrm{stand}})\mathit{\{}\\ {\displaystyle }& {\displaystyle }{S}_{\mathrm{foliage}}={\mathrm{prevS}}_{\mathrm{foliage}}-{\mathrm{EV}}_{\mathrm{foliage}}+R_{\mathrm{tree}\mathrm{or}\mathrm{stand}}\\ {\displaystyle }& {\displaystyle }\mathrm{unintR}=\mathrm{0}\mathit{\}}\\ {\displaystyle }& {\displaystyle }\mathrm{else}\mathit{\{}\\ {\displaystyle }& {\displaystyle }{S}_{\mathrm{foliage}}=C_{\mathrm{foliage}}\\ {\displaystyle }& {\displaystyle }\mathrm{unintR}=R_{\mathrm{tree}\mathrm{or}\mathrm{stand}}-{\mathrm{RemC}}_{\mathrm{foliage}}\mathit{\}}\\ {\displaystyle }& {\displaystyle }\mathrm{else}\mathit{\{}\\ {\displaystyle }& {\displaystyle }{S}_{\mathrm{foliage}}=C_{\mathrm{foliage}}\\ {\displaystyle }& {\displaystyle }\mathrm{unintR}=R_{\mathrm{tree}\mathrm{or}\mathrm{stand}}\mathit{\}}.\end{array}$$

For the leafless period, we have C_foliage=0, which gives unintR=R_{tree or stand}.

Throughfall and stemflow fluxes are then calculated separately for the leaved and leafless periods. For both periods, tree or stand throughfall and pre-stemflow (preSF; in litres) are considered as complementary fractions of the non-intercepted rainfall. Pre-stemflow is the amount of rain deviated towards the branches and the trunk but not necessarily reaching the base of the trunk due to storage and evaporation losses. At the stand level, pre-stemflow is estimated separately for each tree species.

$$\begin{array}{}\text{(47)}& {\displaystyle }& {\displaystyle }{\mathrm{TF}}_{\mathrm{tree}\mathrm{or}\mathrm{stand}}=\mathrm{\%}\mathrm{TF}\cdot \mathrm{unintR}\text{(48)}& {\displaystyle }& {\displaystyle }{\mathrm{preSF}}_{\mathrm{tree}\mathrm{or}\mathrm{sp}}=\mathit{\%}\mathrm{SF}\cdot \mathrm{unintR}\end{array}$$

At this stage, the following state variables are used

• the tree or species bark storage (S_bark; in litres) is the amount of water stored in the bark of a given tree or in that of all the trees of a same-tree species,

• the previous tree or species bark storage (prevS_bark; in litres), which is the tree or species bark storage at the previous time step;

• the remaining bark storage capacity of a given tree or species (RemC_bark; in litres):

$$\begin{array}{}\text{(49)}& {\mathrm{RemC}}_{\mathrm{bark}}=C_{\mathrm{bark}}-({\mathrm{prevS}}_{\mathrm{bark}}-{\mathrm{EV}}_{\mathrm{bark}}).\end{array}$$

Similarly as above for foliage storage and non-intercepted rainfall, various cases are distinguished to update hourly the bark storage and the stemflow volume (SF; in litres) of each tree or species as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{if}({\mathrm{RemC}}_{\mathrm{bark}}>\mathrm{0})\mathit{\{}\\ {\displaystyle }& {\displaystyle }\mathrm{if}({\mathrm{RemC}}_{\mathrm{bark}}>\mathrm{preSF})\mathit{\{}\\ {\displaystyle }& {\displaystyle }{S}_{\mathrm{bark}}={\mathrm{prevS}}_{\mathrm{bark}}-{\mathrm{EV}}_{\mathrm{bark}}+\mathrm{preSF}\\ {\displaystyle }& {\displaystyle }\mathrm{SF}=\mathrm{0}\mathit{\}}\\ {\displaystyle }& {\displaystyle }\mathrm{else}\mathit{\{}\\ {\displaystyle }& {\displaystyle }{S}_{\mathrm{bark}}=C_{\mathrm{bark}}\\ {\displaystyle }& {\displaystyle }\mathrm{SF}=\mathrm{preSF}-{\mathrm{RemC}}_{\mathrm{bark}}\mathit{\}}\\ {\displaystyle }& {\displaystyle }\mathrm{else}\mathit{\{}\\ {\displaystyle }& {\displaystyle }{S}_{\mathrm{bark}}=C_{\mathrm{bark}}\\ {\displaystyle }& {\displaystyle }\mathrm{SF}=\mathrm{preSF}\mathit{\}}.\end{array}$$

At the stand scale, stemflow is obtained by summing stemflow fluxes over the tree species as follows:

$$\begin{array}{}\text{(50)}& {\mathrm{SF}}_{\mathrm{stand}}=\sum _{\mathrm{sp}}{\mathrm{SF}}_{\mathrm{sp}}.\end{array}$$

Tree transpiration

As for evaporation from foliage and bark, the Penman–Monteith equation (see Eq. 27) is used to estimate hourly tree transpiration during the vegetation period. In this case, the radiation absorbed per unit of leaf area by each tree ($\mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}$; in watts per square metre of leaves) is considered and is obtained by

$$\begin{array}{}\text{(51)}& \mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}=\mathit{\%}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}\cdot \mathrm{h}\mathrm{\_}\mathrm{RAD}.\end{array}$$

The individual aerodynamic resistance is determined from Eqs. (36) to (38) applied between the height of largest crown extension (h_lce; in metres) and the dominant height. The individual surface resistance (r_{s_foliage}; in seconds per metre) is defined as the inverse of the foliage stomatal conductance (g_{s_foliage}; in metres per second), which is estimated based on a potential-modifier approach considering soil and climate conditions as well as individual tree characteristics. This approach allows for the increase in stomatal conductance with radiation and for the negative effect of increasing vapour pressure deficit and soil water potential to be accounted for (Granier and Bréda, 1996; Tuzet et al., 2003; Buckley, 2017). For similar soil and climate conditions, the stomatal conductance is acknowledged to be higher for trees with a larger sapwood to leaf area ratio and to decrease with crown height as stomata of top leaves close earlier to avoid cavitation when water stress occurs (Ryan and Yoder, 1997; Schäfer et al., 2000).

$$\begin{array}{}\text{(52)}& r_{\mathrm{s}\mathrm{\_}\mathrm{foliage}}={\displaystyle \frac{\mathrm{1}}{g_{\mathrm{s}\mathrm{\_}\mathrm{foliage}}}}\end{array}$$

$$\begin{array}{}\text{(53)}& \begin{array}{rl}& g_{\mathrm{s}\mathrm{\_}\mathrm{foliage}}=\\ & g_{\mathrm{s}\mathrm{0}\mathrm{\_}\mathrm{foliage}}\cdot {\displaystyle \frac{a_{\mathrm{sapwood}}}{a_{\mathrm{leaf}}}}\cdot {\displaystyle \frac{\mathrm{1}}{h_{\mathrm{lce}}}}\cdot M_{\mathrm{radiation}}\cdot M_{\mathrm{soil}\mathrm{water}}\cdot M_{\mathrm{vpd}},\end{array}\end{array}$$

with g_{s0_foliage} as the reference stomatal conductance (in metres per second), $\frac{a_{\mathrm{sapwood}}}{a_{\mathrm{leaf}}}$ as the sapwood to leaf area ratio (in square metres per square metre) calculated at the tree level (see Jonard et al., 2020a, for details) and M_radiation as the radiation modifier,

$$\begin{array}{}\text{(54)}& {\displaystyle \frac{\mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}}{\mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{tree}\mathrm{\_}\mathrm{leaf}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}+p_{\mathrm{radiation}}}},\end{array}$$

where p_radiation is a parameter characterizing stomatal response to radiation. M_soil water is the soil water modifier,

$$\begin{array}{}\text{(55)}& e^{-p_{{\mathrm{1}}_{\mathrm{SW}}}(p_F-\mathrm{2.5}{)}^{p_{{\mathrm{2}}_{\mathrm{SW}}}}}\mathrm{when}{p}_F>\mathrm{2.5},\mathrm{1}\mathrm{otherwise},\end{array}$$

where p_F (in centimetres) is the base-10 logarithm of the mean soil water potential (ϕ) (mean value of the various horizons weighted based on root proportion; see below in the “Water uptake distribution among soil horizons” section for calculation details of the soil water potential) and $p_{{\mathrm{1}}_{\mathrm{SW}}}$ and $p_{{\mathrm{2}}_{\mathrm{SW}}}$ are two parameters characterizing the stomatal response to soil water potential. M_vpd, the VPD modifier, is

$$\begin{array}{}\text{(56)}& \mathrm{1.0}-p_{\mathrm{VPD}}\cdot \ln \mathrm{VPD},\end{array}$$

where p_VPD is a species-dependent parameter characterizing stomatal response to vapour pressure deficit.

The latent heat flux density (in watts per square metre) determined by applying this parametrization to Eq. (27) is then converted to tree transpiration (TR_tree; in litres per hour) using the same approach as for foliage evaporation that was described in Eqs. (40) and (41). Finally, TR_tree is corrected by multiplying it by the proportion of green leaves (greenProp) and by the fraction of leaves not covered with water ($\mathrm{1}-\frac{S_{\mathrm{foliage}}}{C_{\mathrm{foliage}}})$, considering that transpiration occurs from photosynthetically active and dry leaves only.

Ground vegetation transpiration and soil evaporation

The Penman–Monteith equation is also used to estimate ground vegetation transpiration and soil evaporation at the tree and stand scale. For this purpose, the radiation transmitted to the understory is subdivided for each time step into the radiation absorbed by per unit of leaf area of the ground vegetation ($\mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{grd}\mathrm{\_}\mathrm{veg}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}$; in watts per square metre of leaves) and the radiation absorbed by the soil ($\mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{soil}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}$ in watts per square metre) through application of the Beer–Lambert law:

$$\begin{array}{}\text{(57)}& {\displaystyle }\begin{array}{rl}& \mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{grd}\mathrm{\_}\mathrm{veg}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}=\\ & {\displaystyle \frac{\mathit{\%}\mathrm{transRAD}\cdot \mathrm{rad}\cdot \left(\mathrm{1}-\exp (-k\cdot {\mathrm{LAI}}_{\mathrm{grd}\mathrm{\_}\mathrm{veg}}\cdot {\mathrm{greenProp}}_{\mathrm{stand}})\right)}{{\mathrm{LAI}}_{\mathrm{grd}\mathrm{\_}\mathrm{veg}}\cdot {\mathrm{greenProp}}_{\mathrm{stand}}}}\end{array}\text{(58)}& {\displaystyle }\begin{array}{rl}& \mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{soil}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}=\\ & \mathit{\%}\mathrm{transRAD}\cdot \mathrm{rad}\cdot \exp \left(-k\cdot {\mathrm{LAI}}_{{\mathrm{grd}}_{\mathrm{veg}}}\cdot {\mathrm{greenProp}}_{\mathrm{stand}}\right),\end{array}\end{array}$$

where k is the extinction coefficient fixed to 0.5 (Teh, 2006), LAI_{grd_veg} is the leaf area index of the ground vegetation calculated as the difference between the ecosystem LAI and the tree or stand LAI, and greenProp_stand is the proportion of remaining green leaves at the stand level.

The energy effectively available for soil evaporation is obtained by subtracting the soil heat flux density (G; in watts per square metre) from $\mathrm{h}\mathrm{\_}{\mathrm{aRAD}}_{\mathrm{soil}\mathrm{\_}{\mathrm{m}}^{\mathrm{2}}}$. G is estimated based on the temperature gradient and the soil thermal conductivity (K; fixed to 0.25 W m⁻¹ K⁻¹) as follows:

$$\begin{array}{}\text{(59)}& G=K\times {\displaystyle \frac{T_{\mathrm{surf}}-T_{\mathrm{int}}}{{\mathrm{th}}_{\mathrm{org}}/\mathrm{100}}},\end{array}$$

with T_surf (in degrees Celsius), the temperature at the soil surface, considered as equal to air temperature (T), T_int (in degrees Celsius), the temperature at the interface between the organic layers and the mineral soil (see Jonard et al., 2020a, for more information on the way T_int is obtained), and th_org (in metres), the thickness of the organic layer.

For ground vegetation transpiration and soil evaporation, the aerodynamic resistance is computed by applying Eqs. (36)–(38) between the ground level and the dominant height.

The surface resistances of the ground vegetation (r_{s_grd_veg}) and of the soil (r_{s_soil}) are the reciprocals of the ground vegetation and soil conductances, respectively. The ground vegetation conductance (g_{s_grd_veg}; in metres per second) is estimated based on the same approach as g_{s_foliage} for tree transpiration, while the soil conductance (g_{s_soil}; in metres per second) depends on the relative extractable water (see below for computation details) of the forest floor at the previous time step (prevREW_{forest_floor}) as follows:

$$\begin{array}{}\text{(60)}& \begin{array}{rl}& g_{\mathrm{s}\mathrm{\_}\mathrm{soil}}=\\ & g_{\mathrm{s}\mathrm{\_}\mathrm{soil}\mathrm{\_}\min }+(g_{\mathrm{s}\mathrm{\_}\mathrm{soil}\mathrm{\_}\max }-g_{\mathrm{s}\mathrm{\_}\mathrm{soil}\mathrm{\_}\min })\cdot {\mathrm{prevREW}}_{\mathrm{forest}\mathrm{\_}\mathrm{floor}}.\end{array}\end{array}$$

The latent heat flux density (in watts per square metre) is then converted to ground vegetation transpiration (TR_{grd_veg}; in litres per hour) and soil evaporation (EV_soil; in litres per hour) using the same approach as for tree transpiration and foliage evaporation, Eqs. (40) and (41).

Soil hydraulic properties

The modelling of water uptake distribution among soil horizons and of water transfer from a horizon to another requires estimates of the hydraulic properties for all soil horizons. The relationship between the soil water content (θ; in cubic metres per cubic metre) and the absolute matric potential (h; in centimetres) is described by the van Genuchten function

$$\begin{array}{}\text{(61)}& \mathit{\theta }={\mathit{\theta }}_{\mathrm{r}}+S\cdot ({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}),\end{array}$$

which can be rearranged into the form

$$\begin{array}{}\text{(62)}& S={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}\end{array}$$

and further modified into the form

$$\begin{array}{}\text{(63)}& S={\left[\mathrm{1}+(\mathit{\alpha }\left|h\right|{)}^n\right]}^{-\left(\mathrm{1}-\frac{\mathrm{1}}{n}\right)},\end{array}$$

with θ_r, the residual water content (in cubic metres per cubic metre), θ_s, the saturated water content (in cubic metres per cubic metre), S, the relative water content, and α and n, two parameters.

The Mualem–van Genuchten function allows for the estimation of the soil hydraulic conductivity based on the relative water content and the saturated conductivity.

$$\begin{array}{}\text{(64)}& K=K_{\mathrm{0}}\left(S^{\mathit{\lambda }}{\left\{\mathrm{1}-{\left(\mathrm{1}-S^{n/n-\mathrm{1}}\right)}^{\mathrm{1}-\frac{\mathrm{1}}{n}}\right\}}^{\mathrm{2}}\right),\end{array}$$

with K, the hydraulic conductivity (in centimetres per day), K₀, the saturated conductivity (in centimetres per day), and λ, a parameter.

These two functions (Eqs. 63 and 64) partly share the same parameters, which are estimated based on soil horizon properties (i.e. organic carbon content, C_org, and particle size distribution). For organic horizons, values from Dettmann et al. (2014) are used for α, n and λ (α=0.251; n=1.75; λ=0.5), and the equation of Päivänen (1973) for Sphagnum peat is considered for K₀.

$$\begin{array}{}\text{(65)}& K_{\mathrm{0}}={\mathrm{10}}^{\left(-\mathrm{2.321}-\mathrm{13.22}\cdot {\mathit{\rho }}_{\mathrm{b}}\cdot \frac{\mathrm{1000}}{\mathrm{1}\mathrm{000}\mathrm{000}}\right)\cdot \mathrm{24}\cdot \mathrm{60}\cdot \mathrm{60}},\end{array}$$

with ρ_b the bulk density (in kilograms per cubic metre).

For mineral horizons, pedotransfer equations elaborated by Weynants et al. (2009) are used as follows:

$$\begin{array}{}\text{(66)}& {\displaystyle }\begin{array}{rl}& \ln \mathit{\alpha }=\\ & -\mathrm{4.3003}-\mathrm{0.0097}\cdot \mathrm{clay}+\mathrm{0.0138}\cdot \mathrm{sand}-\mathrm{0.0992}\cdot C_{\mathrm{org}}\end{array}\text{(67)}& {\displaystyle }\begin{array}{rl}& \ln (n-\mathrm{1})=\\ & -\mathrm{1.0846}-\mathrm{0.0236}\cdot \mathrm{clay}-\mathrm{0.0085}\cdot \mathrm{sand}+\mathrm{0.0001}\cdot {\mathrm{sand}}^{\mathrm{2}}\end{array}\text{(68)}& {\displaystyle }\begin{array}{rl}& \ln K_{\mathrm{0}}=\\ & \mathrm{1.9582}+\mathrm{0.0308}\cdot \mathrm{sand}-\mathrm{0.6142}\cdot {\mathit{\rho }}_{\mathrm{b}}-\mathrm{0.1566}\cdot C_{\mathrm{org}}\end{array}\text{(69)}& {\displaystyle }\mathit{\lambda }=-\mathrm{1.8642}-\mathrm{0.1317}\cdot \mathrm{clay}+\mathrm{0.0067}\cdot \mathrm{sand},\end{array}$$

with clay and sand, the clay and sand content of the soil (10⁻² grams per gram), respectively; C_org is the organic carbon content of the soil (in grams per kilogram), and ρ_b is the bulk density (in grams per cubic centimetre).

Water uptake distribution among soil horizons

Once tree and ground vegetation hourly transpiration has been calculated, the module sums the transpiration of all trees for the stand approach and adds the ground vegetation transpiration to obtain the hourly stand transpiration, corresponding to the stand water uptake. Then, tree or stand water uptake is distributed among the horizons according to a method described in Couvreur et al. (2012). This method assumes that water absorption occurs preferentially in horizons where the water potential (matric potential, h, plus a gravimetric component), ϕ, is higher. Moreover, it considers that the amount of water uptake is proportional on the one hand to the difference between the horizon water potential and the averaged water potential weighted by the fine root proportion of the whole soil profile and on the other hand to the fine root proportion of the horizon. This can be transcribed as

$$\begin{array}{}\text{(70)}& \begin{array}{rl}{\mathrm{UP}}_{\mathrm{root}\left(\mathrm{hr}\right)}=& {\mathrm{UP}}_{\mathrm{root}}{f}_{\mathrm{hr}}+K_{\mathrm{comp}}\cdot \mathrm{3600}\left({\mathit{\varphi }}_{\mathrm{hr}}-\sum _{\mathrm{hr}=\mathrm{1}}^N{\mathit{\varphi }}_{\mathrm{hr}}\cdot f_{\mathrm{hr}}\right)\\ & \cdot \mathrm{10}\cdot f_{\mathrm{hr}}\left(A_{\mathrm{pedon}}\mathrm{or}{A}_{\mathrm{stand}}\right)\end{array},\end{array}$$

with UP_root and UP_root(hr),the total water uptake and the water uptake of the “hr” horizon (in litres per hour), respectively, f_hr, the fine root proportion of the horizon hr, K_comp, the compensatory conductivity set to 1.10⁻⁹ (in per second), and ϕ_hr, the horizon water potential (in centimetres).

The second term of the sum of Eq. (70) is null when integrated on all the horizons. Then, it does not change the total amount of water uptake, but it refines its distribution. Moreover, this method can generate water uplift that can occur when the top horizons are much drier than the deep ones.

Water balance of the soil horizons

At the tree and stand scale, the module performs an hourly water balance for each soil horizon hr (numbered from the topsoil) and updates its water content (θ_hr; in cubic metres per cubic metre) as follows:

$$\begin{array}{}\text{(71)}& {\mathit{\theta }}_{\mathrm{hr}}={\mathit{\theta }}_{\mathrm{hr}\mathrm{\_}\mathrm{prev}}+{\displaystyle \frac{({\mathrm{IN}}_{\mathrm{hr}}-{\mathrm{OUT}}_{\mathrm{hr}})}{\mathrm{998}\cdot V_{\mathrm{hr}}}},\end{array}$$

with θ_{hr_prev}, the water content of the hr horizon at the previous time step (in cubic metres per cubic metre), V_hr, the volume of the hr horizon (in cubic metres), IN_hr, the sum of the input water fluxes (in litres), and OUT_hr, the sum of the output water fluxes (in litres).

The input fluxes are the drainage (D; in litres) and the water surplus (S; in litres) from the upper horizon (hr − 1) and the capillary rise (CR; in litres) from the lower horizon (hr + 1) described hereafter and represented in Fig. 1:

$$\begin{array}{}\text{(72)}& {\mathrm{IN}}_{\mathrm{hr}}=D_{\mathrm{hr}-\mathrm{1}}+S_{\mathrm{hr}-\mathrm{1}}+{\mathrm{CR}}_{\mathrm{hr}+\mathrm{1}}.\end{array}$$

The output fluxes are the drainage, the soil evaporation (E_soil; in litres), the root water uptake (UP_root; in litres) and the capillary rise from the current horizon (hr) (Fig. 1).

$$\begin{array}{}\text{(73)}& {\mathrm{OUT}}_{\mathrm{hr}}=D_{\mathrm{hr}}+{\mathrm{EV}}_{\mathrm{soil}\left(\mathrm{hr}\right)}+{\mathrm{UP}}_{\mathrm{root}\left(\mathrm{hr}\right)}+{\mathrm{CR}}_{\mathrm{hr}}\end{array}$$

The water transfer (WT; in litres) between the horizon hr and hr + 1 (considered as drainage if directed downward or as capillary rise if directed upward) is estimated with the Darcy law, and the average conductivity between the horizons is calculated according to the upwind scheme, which takes into account the horizon water potential (e.g. An and Noh, 2014).

$$\begin{array}{}\text{(74)}& \mathrm{WT}={\displaystyle \frac{K_{\mathrm{hr},\mathrm{hr}+\mathrm{1}}}{\mathrm{24}}}\cdot \left({\displaystyle \frac{\mathrm{\Delta }{h}_{\mathrm{m}}}{\mathrm{\Delta }z}}+\mathrm{1}\right)\cdot \left(A_{\mathrm{pedon}}\mathrm{or}{A}_{\mathrm{stand}}\right)\cdot \mathrm{100},\end{array}$$

with

$$\begin{array}{}\text{(75)}& K_{\mathrm{hr},\mathrm{hr}+\mathrm{1}}=\left\{\begin{array}{ll}{K}_{\mathrm{hr}+\mathrm{1}},& {\mathit{\varphi }}_{\mathrm{hr}+\mathrm{1}}>{\mathit{\varphi }}_{\mathrm{hr}}\\ K_{\mathrm{hr}},& {\mathit{\varphi }}_{\mathrm{hr}+\mathrm{1}}\le {\mathit{\varphi }}_{\mathrm{hr}}\end{array}\right.\left(\mathrm{cm}{\mathrm{d}}^{-\mathrm{1}}\right)\end{array}$$

$$\begin{array}{}\text{(76)}& {\displaystyle \frac{\mathrm{\Delta }{h}_{\mathrm{m}}}{\mathrm{\Delta }z}}={\displaystyle \frac{\left|h_{\mathrm{hr}+\mathrm{1}}\right|-\left|h_{\mathrm{hr}}\right|}{\frac{{\mathrm{th}}_{\mathrm{hr}}+{\mathrm{th}}_{\mathrm{hr}+\mathrm{1}}}{\mathrm{2}}\cdot \mathrm{100}}},\end{array}$$

where th (in metres) is the horizon thickness.

To ensure the mass conservation, a variable time step (Δt; in seconds) is considered based on a stability criterion derived from the Peclet number.

$$\begin{array}{}\text{(77)}& \mathrm{\Delta }t={\displaystyle \frac{{\mathit{\theta }}_{{\mathrm{hr}}_{\mathrm{prev}}}\cdot {\mathrm{th}}_{\mathrm{hr}}}{\mathrm{10}\cdot \frac{K_{\mathrm{hr}}}{\mathrm{100}\cdot \mathrm{24}\cdot \mathrm{3600}}}}\end{array}$$

This criterion is calculated for each horizon, and the minimum value is retained. Still, the mass conservation is tested for the whole soil profile at the end of each hour. If the water balance error exceeds 0.01 mm, the time step is divided by 10 (with 1000 as a maximum). The hourly water transfer is then obtained by cumulating the discretized values of water transfer.

For the top horizon, D_hr−1 is initialized at TF + SF, and CR_hr is set to 0. For the current horizon, if WT≥0, D_hr=WT, else D_hr=0 and ${\mathrm{CR}}_{\mathrm{hr}+\mathrm{1}}=-\mathrm{WT}$.

Soil evaporation occurs only in organic horizons. The amount of water evaporated from the horizon hr (EV_soil(hr); in litres) is obtained by taking the minimum value between the remaining water to evaporate (remEV_soil(hr); in litres) and the volume of extractable water in the horizon (${\mathrm{VEW}}_{\mathrm{hr}}={\mathrm{EW}}_{\mathrm{hr}}\cdot \left(A_{\mathrm{pedon}}\mathrm{or}{A}_{\mathrm{stand}}\right)$; in litres). For the upper organic horizon, remEV_soil(hr) is initialized to the total amount of water evaporated from the soil and is progressively decremented by subtractingEV_soil(hr) for the deeper organic horizons as follows:

$$\begin{array}{}\text{(78)}& {\mathrm{remEV}}_{\mathrm{soil}\left(\mathrm{hr}\right)}={\mathrm{remEV}}_{\mathrm{soil}(\mathrm{hr}-\mathrm{1})}-{\mathrm{EV}}_{\mathrm{soil}(\mathrm{hr}-\mathrm{1})}.\end{array}$$

In both mineral and organic horizons, if the water balance leads to a soil horizon water content higher than saturation, the soil horizon water content is set to the value of the saturated water content and a surplus is calculated. Part of this surplus is passed to the next horizon (S_hr−1), while the rest is considered as preferential flows and is added to the deep drainage (DD).

$$\begin{array}{}\text{(79)}& S_{\mathrm{hr}-\mathrm{1}}={\mathrm{IN}}_{\mathrm{hr}}-\left({\mathit{\theta }}_{\mathrm{s}\mathrm{\_}\mathrm{hr}}-{\mathit{\theta }}_{\mathrm{hr}\mathrm{\_}\mathrm{prev}}\right)\cdot V_{\mathrm{hr}}\cdot \mathrm{998}\cdot \left(\mathrm{1}-{\mathit{\nu }}_{\mathrm{hr}}\right)-\mathrm{OUT},\end{array}$$

with v_hr, the additional coarse fraction of the horizon (in cubic metres per cubic metre), not accounted for in the bulk density. The deep drainage is calculated as the sum of D_hr and S_hr−1 of the last horizon plus the preferential flows. Before passing to the next horizon, D_hr−1 takes the value of D_hr and CR_hr the value of CR_hr+1.

Absolute and relative extractable water

The absolute extractable water (EW; in millimetres) is defined as the amount of water stored in the soil that can be used by the plants:

$$\begin{array}{}\text{(80)}& \mathrm{EW}=\sum _{\mathrm{hr}=\mathrm{1}}^n\left({\mathit{\theta }}_{\mathrm{hr}}-{\mathit{\theta }}_{\mathrm{wp}\mathrm{\_}\mathrm{hr}}\right)\cdot {\mathrm{th}}_{\mathrm{hr}}\cdot \left(\mathrm{1}-{\mathit{\nu }}_{\mathrm{hr}}\right),\end{array}$$

where θ_{wp_hr} is the water content of the soil horizon at the wilting point (in cubic metres per cubic metre).

Table 1 Description of the different module parameters for sessile oak and European beech and origin of their value.

The relative extractable water (REW; in millimetres) corresponds to the ratio between this value of extractable water and the reference extractable water at the field capacity (EW_ref; in millimetres):

$$\begin{array}{}\text{(81)}& \mathrm{REW}={\displaystyle \frac{\mathrm{EW}}{{\mathrm{EW}}_{\mathrm{ref}}}},\end{array}$$

with

$$\begin{array}{}\text{(82)}& {\mathrm{EW}}_{\mathrm{ref}}=\sum _{\mathrm{hr}=\mathrm{1}}^n\left({\mathit{\theta }}_{\mathrm{fc}\mathrm{\_}\mathrm{hr}}-{\mathit{\theta }}_{\mathrm{wp}\mathrm{\_}\mathrm{hr}}\right)\cdot {\mathrm{th}}_{\mathrm{hr}}\cdot \left(\mathrm{1}-{\mathit{\nu }}_{\mathrm{hr}}\right),\end{array}$$

where θ_{fc_hr} is the water content of the soil horizon at the field capacity (in cubic metres per cubic metre).

2.2 Parameter determination

Most of the model parameters were taken directly from the literature. In addition, an adjustment of some relationships was conducted using available data, which are described hereafter, but no overall calibration of the model was performed. The model parameters for sessile oak and European beech are presented in Table 1. Regarding common hornbeam, less information is available. For this tree species, we used specific parameters for light interception, photosynthesis, respiration and carbon allocation (Appendix A) but the same parameters as European beech for water balance and phenology given their similar morphology.

For the hydrological module, the parameters of Eq. (53) determining the stomatal conductance were determined based on data from Jonard et al. (2011) using a nonlinear fitting procedure.

For the soil hydraulic properties, the saturation θ_s was based on the 0.999 quantile of measured soil water contents (see Sect. 2.4 for more details). For a horizon without a soil water content sensor, θ_s was interpolated from the closest horizons. Then, the wilting point water content was determined using the obtained saturated water content and the Eq. (63) with a matric potential, h, of 15 000 cm.

The parameters of the phenological module used to calculate the start of budburst were determined using observations from the Pan European Phenology dataset (PEP725), which provides data about phenological observations across different European countries, though not in Belgium. We selected 129 sites on the western border of Germany covering the latitudes of our six study plots (49.5–51.0^∘ N), for which the budburst dates of a representative tree were available at least between 1951 and 2015. The daily minimum, maximum and mean temperatures required to achieve the calibration came from the meteorological stations of the DWD (Deutscher Wetterdienst) Climate Data Center. Phenological data from each site were assigned to the nearest meteorological station (five different stations were sufficient). The calibration was carried out with the Phenological Modeling Platform software (Chuine et al., 2013). This module enables the user to perform a Bayesian calibration procedure using the algorithm of Metropolis et al. (1953). Some of the parameters can also be fixed. In our case, the chilling starting date of the two-phase model was fixed to 1 November of the previous year (e.g. Chiang and Brown, 2007; Roberts et al., 2015) in order to enhance the accuracy of the other parameter calibration. The length of the budburst period (necessary to determine the budburst shift), the leaf development, yellowing and falling rates were all adjusted from phenological observations conducted in our study sites on 20 trees.

2.3 Site description

Six sessile oak (Quercus petraea (Matt.) Liebl.) and European beech (Fagus sylvatica L.) stands located in Wallonia (Belgium) were used to evaluate the model. They all belong to Long-Term Ecological Research sites (Belgian LTER network). Three of them were located in Baileux and were monitored since 2001. The three other stands have been part of the level II plot network of the International Co-operative Programme on Assessment and Forests (ICP Forests) since 1998 and were located in Louvain-la-Neuve, Chimay and Virton. These sites were selected because their contrasted stand structure, species composition, soil and climate make them suitable for testing the ability of the model to account for structure complexity in various ecological conditions (at the regional scale).

2.3.1 Stand characteristics

The experimental site of Baileux was installed to study the impact of species mixture on forest ecosystem functioning (Jonard et al., 2006, 2007, 2008; André et al., 2008a, b) and consisted of three plots. Two plots were located in stands dominated by either sessile oak or European beech, and the third one presents a mixture of both species. In these plots, sessile oak trees originated from a massive regeneration in 1880 and displayed the typical Gaussian distribution of even-aged stands, while European beech trees appeared to be progressively giving rise to an uneven-aged structure with all diameter classes represented. The stand in Chimay was an ancient coppice with standards, presently composed of mature sessile oak trees with an important hornbeam understorey. The stands in Louvain-la-Neuve and Virton were both more or less even-aged stands dominated by European beech but differed in their age, with much older trees in Louvain-la-Neuve than in Virton (130 vs. 60 years old in 2009). All stand characteristics are provided in Table 2.

Table 2 Initial stand characteristics for the main tree species and for the whole stands.

Table 3 Soil and meteorological characteristics of the different study sites (2001–2016 period).

2.3.2 Soil properties

The Baileux, Chimay and Virton stands were all located on Cambisol but with some nuances, ranging from Dystric to the Calcaric variants in Chimay and Virton, respectively, while an Abruptic Luvisol was found in Louvain-la-Neuve (FAO soil taxonomy). All sites presented a moder humus, except Virton for which mull was observed. In Baileux, Chimay and Louvain-la-Neuve, the soil developed from the parent bedrock mixed with aeolian loess deposition that occurred at the interglacial period. In Virton, the soil originated only from the bedrock weathering. The parent materials were sandstone and shales, clayey sandstone, and hard limestone bedrocks in Baileux, Chimay and Virton, respectively. In Louvain-la-Neuve, the soil was almost exclusively built from the loess deposits. These differences in parent material generated contrasting physical and chemical soil properties (Table 3).

The soil textures also varied significantly among sites. Based on the USDA taxonomy, the soil texture was silty clayey loam and silty loam in Baileux and Louvain-la-Neuve, respectively. In Chimay and Virton, finer soil textures were observed, with a clayey loam and a clay texture, respectively. In relation to the texture, drainage was good in Baileux and Louvain-la-Neuve, while the presence of inflating clay triggered the appearance of a shallow water table during the wet period and drought cracks during summer in Chimay. In Virton, despite the high clay content in the lower horizons, drainage was good due to the existence of faults in the bedrock (Table 3).

Finally, stoniness and drainage influenced the estimate of the maximum extractable water reserve. While the beech-dominated and mixed stands in Baileux and in Virton showed the lowest water reserve, the highest value was found in Louvain-la-Neuve, with intermediate values for the oak stand in Baileux and in Chimay (Table 3).

2.3.3 Climate

Even if the same type of climate occurred all over Belgium (temperate oceanic), the study sites were located in different bioclimatic zones (Van der Perre et al., 2015). Louvain-la-Neuve was in the Hesbino-Brabançon zone with the highest average temperatures (11.0 ^∘C) between 2001 and 2016 and the driest conditions (818 mm). Despite their close locations, Baileux and Chimay were part of different zones. Baileux was in Basse et moyenne Ardenne, while Chimay was in the zone Fagne, Famenne et Calestienne. Average temperatures are similar for both locations (i.e. 9.8 ^∘C in Baileux and 9.7 ^∘C in Chimay). Yet, a consistent difference in terms of precipitation is observed. Baileux being more elevated, it receives on average 1075 mm of precipitation each year, while only 940 mm is measured in Chimay with respect to the rainfall–altitude relationship (Poncelet, 1956). Finally, Virton was part of the Basse Lorraine with elevated annual rainfall (1060 mm) and intermediate average temperature values (9.9 ^∘C) (Table 3).

For Chimay, Louvain-la-Neuve and Virton, we used data from the meteorological stations of the Pameseb network. The records covered the 1999–2018 period. A tipping bucket located at 1 m height was used to monitor rainfall. Global radiation was registered with a pyranometer; air temperature was registered with a resistance sensor thermistor; relative humidity was registered with a psychrometer; wind speed was registered with an anemometer. All of these devices were placed at 1.5 m height. Data were collected at 12 min intervals and were then averaged hourly. For Baileux, an independent meteorological station managed by our laboratory has been used to collect meteorological data since 2002. The devices were identical to those described before. Air temperature, relative humidity and rainfall were monitored at 1.5 m. Wind speed and global radiation were taken at 2.5 m above the ground.

2.4 Model evaluation

The various routines to calculate the budburst starting date were tested, and the two-phase model with the optimum response function for chilling was retained for the evaluation of the water balance module, as this approach performed better (see Sect. 3.1.1).

2.4.1 Phenology

The phenological observations available on the level II sites of Chimay, Louvain-la-Neuve and Virton were used to evaluate the model predictions. These phenological observations were carried out on 20 dominant and co-dominant sessile oaks in Chimay (2012–2014) and 20 dominant and co-dominant European beeches in Louvain-la-Neuve and Virton (2012–2016) according to the ICP Forests manual (Beuker et al., 2016). They consisted of weekly observations of the percentage of budburst, yellowing and leaf fall, depending on the season. As the model predicted the budburst for an average tree, we evaluated it with the budburst observations of the median tree. In addition, we visually assessed the agreement between the predicted and observed increase in leaf biomass proportion (leafProp) during the leaf development period and between the predicted and observed decrease in green leaf proportion (greenProp) and in leafProp during leaf yellowing and leaf fall, respectively. We did not perform a statistical evaluation for these latter variables as the corresponding processes were not calibrated independently in the model. Finally, as there were no data available for trees of different social status, we could not directly evaluate the option Phenology at tree level. We evaluated however its impact on tree growth predictions for the three stands in Baileux.

2.4.2 Water balance

Regarding the water balance module, the evaluation was conducted using variables integrating most of the processes described in the model. The observed throughfall, extractable water dynamics, individual transpiration and deep drainage (considered in the next section) were compared to model predictions. For the evaluation of the throughfall, extractable water and drainage predictions, we used simulations carried out at the stand scale since the corresponding observations cannot be related to a particular tree. Regarding individual tree transpiration, the approaches at the two scales were compared (tree vs. stand).

For the evaluation of throughfall predictions, only independent throughfall data collected in Chimay, Louvain-la-Neuve and Virton between 2000 and 2016 were used, as the rainfall partitioning routine was calibrated based on data from the Baileux forest (André et al., 2008a, b). The collecting devices consisted of three long gutters connected to plastic barrels. The throughfall volume was measured weekly based on the height of water in the barrels. A log transformation of both the observations and the predictions was necessary to remove the heteroscedasticity.

Individual tree transpiration predictions were evaluated against observations derived from sap flux measurements. These measurements were taken on 16 sessile oak and 16 European beech trees of different sizes in the three stands of Baileux between April and September 2003 (Jonard et al., 2011).

Extractable water was estimated based on Eq. (80) using soil water content measurements taken between 2005 and 2017 in Baileux and for the 2015–2018 period in the other sites. Soil water content was measured hourly using a thermo-time domain reflectometer (TDR) inserted in some horizons (measurements at three to five different soil depths depending on the site). In order to decrease the influence of the soil disturbance due to the instrument installation, the first year of records was discarded. Indeed, Walker et al. (2004) showed that inserting a moisture sensor in a soil disturbed its hydraulic properties and water content during at least 9 months. The electrical signal from the TDR was transformed into relative dielectric permittivity and then converted into soil volumetric water content (in cubic metres per cubic metre) using the equation of Topp et al. (1980) for Baileux and resorting to our own calibration for the other sites (established based on gravimetric measurements of soil water content).

2.4.3 Drainage

Deep drainage can represent a large water output but is difficult to measure directly. Among the existing indirect approaches to estimate this component, we retained the mass-balance method using chloride ion (Cl⁻) as a tracer. This method has been widely used to estimate groundwater recharge (e.g. Bazuhair and Wood, 1996; Ting et al., 1998; Scanlon et al., 2002) but can be applied to assess deep percolation as well (Willis et al., 1997). It relies on the fact that Cl⁻ is not subject to any chemical transformations in the soil and undergoes only temporary storage in soil (Öberg, 2003). The only Cl⁻ input in our study plots comes from throughfall and stemflow and can be determined from Cl⁻ deposition data obtained from monthly chemical analyses of throughfall and stemflow samples. For the deep drainage, which constitutes the only output, the Cl⁻ concentration is also obtained from monthly chemical analyses of soil solution collected with zero-tension lysimeters at 1 m depth in the three stands of Baileux between 2008 and 2016 and between 2013 and 2016 for the other sites. Deep drainage was estimated yearly by considering that the Cl⁻ amount leaving the soil through drainage was equal to the Cl⁻ input from throughfall and stemflow. As there is a clear annual pattern with a recharge and a discharge period in our study sites, the annual time step is therefore required to verify the hypothesis that chloride concentration in rainfall and in the soil are in a steady-state balance. Based on Eq. (83), the deep drainage flux was estimated and compared to our predictions.

$$\begin{array}{}\text{(83)}& \mathrm{Drainage}=(\mathrm{Throughfall}+\mathrm{Stemflow})\cdot {\displaystyle \frac{[\mathrm{Cl}{]}_{\mathrm{Througfall}-\mathrm{Stemflow}}}{[\mathrm{Cl}{]}_{\mathrm{Drainage}}}},\end{array}$$

with [Cl]_{Throughfall−Stemflow}, Cl⁻ concentration in throughfall and stemflow, and [Cl]_Drainage, Cl⁻ concentration in drainage water.

2.4.4 Statistical analyses

To test the quality of the predictions, different statistical tests and indices were used. The absolute bias, defined as the difference between the mean observation and prediction, and the relative bias, corresponding to the ratio between the absolute bias and the mean observation, were calculated to detect any over- or underestimation. To assess the precision of the predictions, the root mean square error (RMSE) was used and calculated as follows:

$$\begin{array}{}\text{(84)}& \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\sum {\left(\mathrm{obs}-\mathrm{pred}\right)}^{\mathrm{2}}}{n}}},\end{array}$$

with n the number of observations.

When the range of values differed considerably for one variable between the different sites, the RMSE was divided by the range, i.e. the difference between the maximum and the minimum values. This normalized root-mean-square error (NRMSE) is much more adapted for comparisons in these situations.

The agreement between observations and predictions was also evaluated with the Pearson's correlation coefficient (r) and with a regression test conducted to analyse the linear relationship between observed and predicted values. As both predictions and observations are subject to uncertainties, we used Deming regression. Then, we tested whether the regression line confidence interval (95 %) included the identity line. These tests were realized with the mcr package in R.

3 Results

3.1 Evaluation of model performance

3.1.1 Phenology

On average, the budburst was best predicted with the two-phase model with the optimum response function for chilling (bias = 2.46 d compared with 8.23 and −5.88 d for the one-phase and two-phase with sigmoid response function for chilling models, respectively). However, this option was less appropriate to capture the interannual variations (Pearson's r=0.537) than the one-phase model (Pearson's r=0.680). The temporal variability was very poorly estimated with the two-phase model with the sigmoid response function for chilling, which displayed an inverse trend for the ranking among years (Pearson's $r=-\mathrm{0.277}$) (Fig. 2a). Moreover, as this latter model was not able to predict the end of the chilling period for some years in Louvain-la-Neuve (European beech), all results for this site were discarded. The predicted leaf development displayed a good agreement with observations (Fig. 3).

Figure 2 Comparison of the observed and predicted budburst of the median tree in Chimay, Virton and Louvain-la-Neuve for the three phenological variants implemented. The quality of the predictions is indicated by the RMSE, the absolute bias and the Pearson's correlation coefficient (r).

Figure 3 Observed and predicted increase in leaf proportion in Chimay, Louvain-la-Neuve and Virton during the budburst and leaf development phase (data from 2012 to 2016). Observations are missing in Chimay for 2013, in Louvain-la-Neuve for 2012 and 2013, and in Virton for 2013.

Simulated leaf yellowing and leaf fall were also evaluated by comparison with observations. While the leaf ageing threshold was taken from Dufrêne et al. (2005), the yellowing parameter determining the length of the yellowing period was adjusted with the 5 years of data from Chimay (sessile oak), Louvain-la-Neuve and Virton (European beech). Therefore, only the yellowing start was independently evaluated. The prediction of the start of the yellowing displayed a low absolute bias (2.7 d) and RMSE (7.0 d). However, a weak correlation (0.056) was found between predictions and observations (data not shown).

For the temporal dynamics of leaf yellowing and leaf fall, the agreement between model predictions and observations was just assessed visually since the parameter regulating these processes (yellowing, falling rate and falling frost amplifier) were adjusted with the same data. The overall agreement was good. The simulated decrease in green leaf proportion was similar for all sites as the photoperiod reduction is identical for each site and year (Fig. 4a, c and e). The only noticeable difference came from the yellowing starting date, which depended on air temperature. For Chimay (sessile oak), a close agreement was found between predictions and observations. For Louvain-la-Neuve (European beech), predictions were correctly centred, but the predicted trend was more abrupt, and the start of the decrease displayed some delay, except in 2012. For Virton (European beech), the decreasing trend was correctly displayed, but the decrease start was less precise in 2016 (Fig. 4e).

Figure 4 Observed and predicted temporal dynamics in leaf yellowing and in leaf fall in Chimay, Louvain-la-Neuve and Virton (data from 2012–2016). Yellowing is represented by the decrease in green leaf proportion (a, c, e) and leaf fall by the decrease in total leaf proportion (b, d, f).

Concerning the leaf fall, the temporal dynamics were correctly represented in Chimay (sessile oak). In Louvain-la-Neuve (European beech), the model predicted a slightly too slow decrease in leaf proportion in 2012 and 2015. For the other years, the observed and predicted leaf proportion matched well even if the predicted start of the fall appeared later than in the observations for some years. In Virton (European beech), the predictions were well centred with regard to the observations, but the decrease in leaf proportion was a bit too fast in 2012 (Fig. 4b, d and f).

The option phenology at tree level was used to test if the agreement between predicted and observed basal area increment could be improved. With the default phenology option, HETEROFOR tended to overestimate the growth of dominant trees and underestimate that of suppressed trees (Jonard et al., 2020a). With the option phenology at tree level, this bias was partially resorbed. The slope of the Deming regression went from 0.74 to 0.84 for sessile oak and from 0.79 to 0.88 for European beech, being much closer to the identity line (Appendix D). This was however at the price of slightly lower Pearson's correlation coefficients.

3.1.2 Water balance

For each site, the main water fluxes affecting the water balance were calculated daily, summed up and the annual values were averaged for the 2002–2016 period (Table 4). Depending on the site, 65 % to 78 % of the rainfall reached the floor as throughfall and 6 % to 13 % as stemflow. The remaining 16 % to 22 % was intercepted by the tree foliage and the bark and evaporated. Then, 31 % to 45 % of the water received as rainfall returned in the atmosphere through tree transpiration. The remaining 26 % to 44 % was lost from the ecosystem through drainage.

Figure 5 (a) Comparison of the log-transformed observed and predicted monthly throughfall in Chimay (oak), Louvain-la-Neuve (beech) and Virton (beech) between 2000 and 2016 and comparison of observed and (b) predicted annual drainage in all study stands between 2008 and 2016. The quality of the non-transformed predictions is indicated by the RMSE, the relative bias and the Pearson's correlation coefficient (r). The shaded area represents the confidence interval of the Deming regression (95 %) of observations on predictions, and the solid line corresponds to the identity line.

3.1.3 Rainfall partitioning

Rainfall partitioning was correctly reproduced by the HETEROFOR model. Across all considered sites (Virton, Chimay and Louvain-la-Neuve), the mean bias of throughfall predictions was very limited (−1.3 %) and nonsignificant (P value of the paired t test is 0.316). The confidence interval of the linear relationship between the logarithm of the observed and predicted throughfall contained the identity line corresponding to the perfect match (upper part of Fig. 5). The correlation between predictions and observations amounted to 0.86, and the RMSE was 16.62 mm, which corresponded to 34.2 % of the mean throughfall (48.6 mm). The separate examination of the different sites revealed that throughfall in Virton (European beech) were very well predicted but that a slight underestimation of the throughfall predictions in Chimay (sessile oak) was compensating for an overestimation of similar magnitude in Louvain-la-Neuve (European beech) (Appendix B).

Figure 6 Comparison of the observed and predicted daily transpiration of sessile oak and European beech in 2003, considering the tree and the stand scale for the water balance calculation. The quality of predictions is indicated by the RMSE, the relative bias and the Pearson's correlation coefficient (r). The shaded area represents the confidence interval of the Deming regression (95 %) of observations on predictions, and the solid line corresponds to the identity line.

Figure 7 Temporal dynamics of observed and predicted extractable water amount (in millimetres) in the various stands. The prediction quality is indicated by the NRMSE, the relative bias and the Pearson's correlation coefficient (r).

3.1.4 Transpiration

The model reproduced well the individual transpiration for sessile oak and European beech in the Baileux site (in 2003), with similar performances at the tree and stand scale (Pearson's r of 0.84–0.85 for sessile oak and 0.88–0.89 for European beech). For European beech, the tree approach corrected the slight bias observed with the stand approach due to an overestimation of high transpiration values. Regarding sessile oak, the small underestimation of transpiration remained whatever the scale considered (Fig. 6).

3.1.5 Soil water content

As the temporal variation in the extractable water was affected by all the water fluxes, it was used to check the performances of the water balance module (Fig. 7). A clear seasonal pattern appeared. At the beginning of the vegetation period, the extractable water (EW) values were highest. Then, tree and ground vegetation transpiration progressively depleted the water reserve, which was partly recharged with rainfall events. Depending on their frequency, duration and intensity, the decline in EW was more or less pronounced, and available water could reach levels close to zero. For all the sites, the Pearson's correlation between observed and predicted relative extractable water ranged from 0.893 to 0.950. These high correlation values and the graph inspection show that the seasonal pattern was precisely reproduced by the HETEROFOR model. NRMSE values range from 10.54 to 13.96 %, while relative bias values were around −2 % and −3 % in the sites Baileux-oak, Baileux-mixed and Chimay and close to −8 % and −9 % in Baileux-beech, Louvain-la-Neuve and Virton. These higher negative biases in the latter stands originated mainly from the model underestimation of the high values of EW (i.e. during wet periods). Despite these similar statistical results, the amount of extractable water in Virton displayed some peculiarities with regard to the other stands. Indeed, the observed EW levels fluctuated considerably more than in the other sites, with frequent peaks both for high and low values that were not represented by the model. Finally, apart from Virton where some discrepancy between observations and predictions can be pointed out, the model quality did not decrease in Chimay or Virton during the 2018 summer, which was categorized as exceptionally dry by the Royal Meteorological Institute of Belgium. The comparison of the tree and stand approach in 2003 indicates that the extractable water calculated at the tree scale progressively deviated from that obtained at the stand scale during the course of the vegetation period and became slightly lower, especially in the sites Baileux-oak and Baileux-mixed (Appendix C). On these graphs, one may notice the heterogeneity in extractable water within the various stands.

3.1.6 Drainage

The predicted deep drainage was compared with estimates calculated on a yearly basis using Cl as a tracer. The RMSE (100.6 mm) and the bias (−19.9 %) were quite large, but a surprisingly good correlation was found between the predicted and estimated drainage (Pearson's r=0.963). Due to the systematic bias, the identity line was not within the confidence interval of the Deming regression despite a regression slope of 0.97 (lower part of Fig. 5).

4 Discussion

In order to predict the impact of global changes on forests, it is crucial to integrate and structure the existing knowledge in process-based models. However, this first step is not sufficient. A detailed documentation of the models and an evaluation of their performances are also needed in order to use them while knowing exactly their strengths and limits. While most models were described in scientific articles or reports, their evaluation was often limited to one or two sites used to illustrate the model functioning and were generally based on integrative response variables such as radial tree growth (Vanclay and Skovsgaard, 1997; Schmidt et al., 2006). Yet, to provide robust predictions of tree growth under changing conditions, the model must be able to accurately reproduce not only the observed tree growth but also the intermediate processes describing resource availability (light, water and nutrients) (Soares et al., 1995). In the following section, we discuss the quality of the predictions for two main drivers of tree growth (phenology and water balance) in relation with the concepts used to describe them.

4.1 Phenology

The two-phase model with the optimum response function for chilling was the least-biased variant for predicting budburst. However, the one-phase model including only the forcing period better captured the interannual variability. While the bias is likely to originate from the model calibration (data used for calibration were observations from western Germany) and could be corrected, the ability of the model to predict temporal variability is more representative of its structural quality. It is common that models accounting only for the forcing period better represent the observations of budburst temporal variability (Leinonen and Kramer, 2002; Fu et al., 2014; Basler, 2016). Two reasons can explain this. First, in areas where the chilling requirements are always met as in western Europe, the chilling parameters are not constrained enough by the observations and therefore are more difficult to calibrate. When few phenological observations are available, taking the chilling into account increases the model complexity without improving the accuracy of the predictions (Leinonen and Kramer, 2002; Fu et al., 2014). Secondly, two-phase models predict the break of endodormancy (end of chilling) and the start of budburst, while they are only calibrated from budburst observations. When endodormancy break observations are used for the calibration, the chilling parameters are estimated more accurately (Chuine et al., 2016), but these data are difficult to obtain and consequently very scarce (Jones et al., 2013).

Similarly to Chuine et al. (2016), we advise using a two-phase model (with the optimum response function for chilling) when endodormancy break observations are available or for long-term simulations with climate conditions beyond the range used for the model calibration. For example, for trees located at the southern margins of their species distribution area, the expected rise in winter temperature could prevent the fulfilling of their chilling requirements and inhibit budburst (Clark et al., 2014). This phenomenon would only be captured by a two-phase model. When few phenological data are available for calibration, the one-phase model should be preferred, especially for simulations with climate conditions similar to those used for the calibration. This highlights the importance of having both one and two-phase options to describe budburst. Most process-based models listed in Table 5 had however only one phenological variant except 4C.

Table 4 Predicted annual water fluxes and the corresponding percentage of rainfall in brackets for the different study sites during the period 2002–2016. The minimum, maximum and mean values from literature are indicated with the number of studies (n) they are based on. The studies taken into account were restricted to sites dominated by beech or by oak in temperate regions with similar meteorological conditions. Data from the same site were averaged so that long monitoring studies do not influence the average value too much.

Papers included in the literature review: Cepel (1967), Aussenac (1968, 1970), Lemée (1974); Nagy (1974); Szabo (1975); Aussenac and Boulangeat (1980); Matzner and Ulrich (1981); Rowe (1983); Bücking and Krebs (1986); Gerke (1987); Giacomin and Trucchin (1992); Neal et al. (1993); Leuschner (1994); Ulrich et al. (1995); Heil (1996); Tarazona et al. (1996); Bellot and Escarre (1998); Didon-Lescot (1998); Herbst et al. (1998); Nizinski and Saugier (1998); Forgeard et al. (1980); Granier et al. (2000); Bent (2001); Michopoulos et al. (2001); Knoche et al. (2002); Mosello et al. (2002); Dripps (2003); Bastrup-Birk and Gundersen (2004); Hanson et al. (2004); Ladekarl et al. (2005); Schipka et al. (2005); Vincke et al. (2005); Carlyle-Moses and Price (2006); Christiansen et al. (2006); Roberts and Rosier (2006); Schmidt (2007); Herbst et al. (2008); Staelens et al. (2008); Ahmadi et al. (2009); Müller and Bolte (2009); Risser et al. (2009); and Gebauer et al. (2012).

A possible improvement of the two-phase phenological models would be to integrate the photoperiod effect on budburst. Indeed, some recent studies have shown evidence that photoperiod can compensate for a lack of chilling temperature that would prevent the buds from opening and for an early frost episode that would trigger budburst before winter (Vitasse and Basler, 2013; Pletsers et al., 2015). This mechanism is particularly present for late-successional species like beech and oak trees and is regularly cited as a key element used to simulate the phenology under climate change (Basler and Körner, 2012). Some models tried to account for the photoperiod effect simply by replacing chilling by photoperiod (Kramer, 1994; Schaber and Badeck, 2003) but, in this way, failed to represent the combined effect of these variables. Recently, a few models integrating the compensatory effect of photoperiod on chilling have appeared. However, these models include more phenological parameters for similar predictive ability (Gauzere et al., 2017). It remains indeed difficult to disentangle the co-varying effect of chilling and day length with in situ measurements (Flynn and Wolkovich, 2018) since photoperiod variations only occur for sites with different latitudes where other confusing factors play a role as well (Primack et al., 2009). Therefore, a large quantity of data are necessary to calibrate these models. Therefore, we decided to favour the accuracy of our phenological model with a more process-based approach, but we are looking forward for improvements in these kinds of models and a more consensual body of literature.

The better growth predictions obtained for the small trees when the phenology was calculated at the individual scale highlight the importance of the “phenological avoidance strategy” displayed by understory trees. This had already been mentioned by Lopez et al. (2008) who observed that early-leafing species received between 45 % and 80 % of their photon flux during the budburst period. Moreover, a simulation study showed that a 1-week (2-week) lengthening of the understory vegetation period with regard to the overstorey in both spring and autumn generated a productivity increase of 32 % (55 %) on such a short period (Jolly et al., 2004).

4.2 Water balance

In a first step, the annual water fluxes predicted by HETEROFOR were compared to measurements and predictions of other studies (Table 4). Then, some water fluxes were individually evaluated when data were available. Finally, some potential improvement of the water balance module were discussed.

Various studies were taken from the literature to compare our water module predictions with observations. They cover a range of annual rainfall comprised of between 425 and 1476 mm (Table 4), which is comparable to what can be found in Belgium. The proportions of rainfall converted to stemflow obtained with HETEROFOR (6.1 % to 13.1 %) are within the range reported in the literature (0.6 % to 20.4 %). This large observation spectrum comes from the important seasonal (higher stemflow proportion in winter than during the vegetation period) and species differences (stemflow importance is higher for beech than oak trees), which are accounted for in HETEROFOR. However, the mean value from the literature (7.3 % of rainfall) is close to the average value for the six study sites (10.3 %). The proportions of intercepted rainfall (15.9 % to 22.0 %) and throughfall (64.8 % to 78.0 %) are also consistent with the ranges reported in other studies (1.9 % to 31.0 % and 59.8 % to 83.1 %). Moreover, we observed a good matching between the average values (respectively 19.5 % and 73.8 % from literature and 19.4 % and 70.2 % for our study sites). For transpiration, the range found in the literature is large (14.8 % to 52.3 %, with an average value of 31.9 %), which is not surprising since interannual and intersite variabilities are high for this variable (Schipka et al., 2005; Vincke et al., 2005). The predicted transpiration proportions are less variable (31.2 % to 44.9 %) and their average value of 36.0 % is slightly above to the mean observed transpiration (31.9 %). Regarding drainage, no direct measurements can be made; all the estimates from the literature come from indirect methods or modelling also subject to uncertainties. The range of drainage values reported in the literature (13 % to 70 %) is very large and contains the range obtained with HETEROFOR (26.3 % to 44.2 %). The mean predicted drainage (39.7 %) is close to the mean value of the literature (37.5 %). By this comparison with the water fluxes reported in the literature, we show that HETEROFOR provides plausible estimates of the various components of the water cycle.

Comparing predicted and observed throughfall is interesting for the evaluation of the water balance module since throughfall is an integrative variable depending on the water storage capacity of foliage, on evaporation and on the proportion of stemflow. The good agreement between observations and predictions indicates that the partitioning of rainfall when passing through the canopy and the evaporation of the water intercepted by foliage and bark are described well. Among the different models of Table 5, no one accounts separately for stemflow and throughfall, but other models not included in the list consider the two fluxes separately (e.g. Gotilwa+ and Castanea). Yet, separating throughfall and stemflow is important, especially for structurally complex stands. In these stands, rainfall interception cannot be simulated based on a mean foliage storage capacity and a mean partitioning between throughfall and stemflow since these parameters vary with stand composition and structure. Our tree-level approach estimating foliage storage capacity and stemflow proportion based on individual tree characteristics allows us to overcome this difficulty. Moreover, if one wants to accurately describe the nutrient cycle, partitioning rainfall is essential as nutrient concentrations in stemflow and throughfall can be 10 to 100 times higher than in rainfall due to dry and wet deposition and canopy exchange (Levia Jr. and Herwitz, 2000; André et al., 2008c; van Stan and Gordon, 2018). Even if the rainfall partitioning can still be improved from a theoretical perspective (e.g. including canopy drainage after rain events or the impact of wind on the foliage storage capacity like in Hörmann et al., 1996, or Muzylo et al., 2009), we chose to limit the level of complexity in order to avoid calibration difficulties.

HETEROFOR satisfyingly reproduced individual tree transpiration with similar prediction quality for the tree and stand approach regarding the water balance calculation. For European beech, the water balance calculation at the tree scale allowed even the correction of the small bias which appeared with the stand approach (Fig. 6). The year selected for the simulation (2003) was particularly dry and hot during summer, which allowed a large range of meteorological and soil water conditions to be covered. It is indeed interesting to test the tree approach under dry conditions since horizontal water redistribution is much less efficient in this case.

Twenty to thirty percent of the transpiration variability remained unexplained by the model, which can be partly ascribed to model inaccuracies but also to the large uncertainty associated with the sap flux measurements. Among other things, the measurements made by Jonard et al. (2011) did not take the azimuthal variation in the sap flux into account since only one sensor per tree was installed.

This first evaluation of tree transpiration predictions indicates that no loss of precision occurred with the tree scale approach, while this detailed spatial representation could have increased the variability in transpiration predictions since it generated some heterogeneity in soil water availability (Appendix C). These good results show that the water balance calculation at the tree scale provides a promising tool to better understand the individual variability and local environment effects on tree water use and sensitivity to drought. This must be considered in a dynamics of continuous improvement of the model and will require more transpiration measurements and in-depth comparisons of predictions and observations.

The amount of extractable water (EW), directly influenced by tree transpiration and soil evaporation, is also a key element of the water cycle, driving, among other things, the drought resistance of a stand. The temporal dynamics of EW was captured well by HETEROFOR, as evidenced by the high correlations (Pearson's coefficient was between 0.893 and 0.950) between observed and predicted EW for the various study sites (Fig. 7). These correlations are within the high end of the range reported for similar models. With the BALANCE model, Grote and Pretzsch (2002) obtained a Pearson's correlation of 0.85 between the observed and predicted soil water content of the upper soil (0–20 cm horizon) in a beech forest in Germany (Freising). Applying BALANCE on three broadleaved stands of oak or beech in Germany, Rötzer et al. (2005) were also able to correctly reproduce soil water content dynamics but mentioned a significant decrease in the quality of predictions during the 2003 drought due to an overestimation of the soil drying, which was not observed with HETEROFOR in 2018. Comparing the observed soil water content at various soil depths with that predicted by the 4C model in mixed oak and pine forest (Brandenburg, Germany), Gutsch et al. (2015) obtained Pearson's correlations ranging from 0.59 to 0.74. In an oak stand in Tennessee (USA), Hanson et al. (2004) compared the ability of nine process-based forest models to reproduce soil water dynamics in the 0–35 cm horizon of the soil and obtained correlations ranging from 0.81 to 0.96.

In the study of Hanson et al. (2004), relative bias was evaluated as well for soil water content and ranged between −1.3 % and 4.0 %. These values are comparable to those found in this study yet a bit lower. Furthermore, discrepancies between predicted and observed EW occurred during limited periods. Several reasons can be advanced to explain them. Errors in the prediction of the budburst date can result in a too early or too late restarting of tree transpiration and induce an inaccurate depletion of the soil extractable water during the vegetation period. In order to distinguish this error source from the others, one could force the model with the observed budburst date. This option is however not yet implemented in the model. The lack of agreement between observed and predicted EW could also be ascribed to the strong heterogeneity of soil properties in forest ecosystems. Similarly, local rainfall events recharging soil extractable water during summer (often associated with thunderstorms) are sometimes not correctly taken into account when missing meteorological data (due to failed sensors or other technical problems) are replaced by rainfall data of a meteorological station further away.

Table 5 Comparison of the spatial scale (S: stand; C: cohort; I: individual; and I^*: individual target tree) and concepts used for describing phenological and hydrological processes in HETEROFOR and in other individual and cohort-based models. Backslash is used to distinguish the various model options. Abbreviations used in for describing transpiration (P–M: Penman–Monteith; SPAC: soil–plant–atmosphere continuum).

^a Grote and Pretzsch (2002); ^b Rötzer et al. (2010); ^c Friend et al. (1997); ^d Seidl et al. (2012); ^e Duursma and Medlyn (2012); ^f Duursma (2008); ^g Simioni et al. (2016); ^h Gutsch et al. (2015); ⁱ Model description on 4C website; ^j Deckmyn et al. (2008); ^k Grote et al. (2011); ^l Collalti et al. (2014); ^m Collalti et al. (2016).

Simplifications and errors in the model conception may also generate divergence between observations and predictions. However, this structural uncertainty can be limited by selecting the most appropriate concepts. HETEROFOR predicts water transfer between soil horizons using the Darcy law. We tried to implement an approach of intermediate complexity between simple bucket models and the Richards equations. From a theoretical point of view, the Richards approach is the most state-of-the-art but requires very long calculation times (Fatichi et al., 2016) and is usually implemented in models specifically dedicated to water flow simulations (in Table 5, only one of the models, MAESPA, uses them). Forest ecosystem models generally use simpler approaches such as the bucket model declined in a large variety of forms (Table 5). These models consider one or several buckets with a specified water storage capacity that is filled with rainfall and is emptied by evapotranspiration. If the soil water content is at field capacity, water is transferred to the underlying layer and finally lost by drainage. Improved versions can account for transfer between buckets in unsaturated conditions using the Darcy law (leaky bucket model).

Our water transfer routine discretizes the soil in horizons whose thickness varies from a few centimetres (upper horizons) to half a metre (deeper horizons). Compared to the numerical resolution of Richards equation, which requires thin soil layers (1 to 2 cm), our vertical discretization of the soil profile is quite coarse and inaccurately predicts the advance of the wetting front. As the tree transpiration and photosynthesis depend on the soil water conditions of the whole soil profile, this inaccuracy has very limited implications on the simulated tree growth. In our approach, water transfer during a time step is calculated based on the horizon water potentials estimated at the end of the previous time step. As such, the model makes the hypothesis that the water content does not change significantly during the time step, which is certainly not the case close to the wetting front and cannot ensure mass conservation. In order to limit this problem, the model uses an adaptive time step estimated based on the Peclet number described in Eq. (77). This allows mass conservation to be ensured.

Finally, another reason that could explain the discrepancy between predictions and observations is the presence of macropores that cause preferential flows. These water fluxes defined as water movements in the soil along preferred pathways that bypass the soil matrix (Hardie et al., 2011) can be generated by soil shrinkage, root growth, chemical weathering, cycles of freezing and thawing or bioturbation (Aubertin, 1971). These macropores are more frequent in forest soils than in agricultural soils as the latter are often ploughed and homogenized. They are however difficult to characterize given their strong spatial heterogeneity in both vertical and horizontal directions (Aubertin, 1971). Adaptations of the Richards equations can be used to account for the preferential flows (dual porosity and dual permeability) but require a good characterization of soil macropores (not possible to achieve routinely in forest soils given their heterogeneity) and are still more complicated to solve than the classical Richards equations. We implemented in the model the transfer of the soil water surplus (when water saturation is reached) to the underlying horizon and the possibility to redirect part of this surplus as deep drainage to account empirically for preferential flows. Indeed, preferential flows in macropores become significant only when rainfall exceeds the water infiltration rate in the soil matrix and accumulates in the soil surface. The fraction of the water surplus considered as preferential flows is an empirical parameter reflecting the macroporosity of the site.

The performances of the soil water transfer routine can also be checked based on the deep drainage flux. In this study, we compared the deep drainage estimated with HETEROFOR and with the chloride mass balance approach. The mean drainage predicted with HETEROFOR was 379 mm yr⁻¹ while the average drainage obtained with the chloride approach amounted to 472 mm yr⁻¹, which corresponds to a bias of −19.9 %. The correlation between the two types of estimate amounted to a Pearson's coefficient of 0.963, with a RMSE value of 100.6 mm. These values depict a constant negative bias in the predictions that can easily be seen on the lower part of Fig. 5. It is hard to tell whether the gap originates from the model or the method used to estimate drainage from the chloride approach. It is more likely that the bias must be ascribed to both. Indeed, on the one hand, even if the use of chemical tracers to estimate drainage or groundwater recharge is commonly used (Scanlon et al., 2002), its application remains subject to uncertainties. First, the chloride method supposes that the main chloride source is rainfall and that the other sources can be neglected (Murphy et al., 1996). This hypothesis is not always fulfilled due to anthropogenic chloride introduction (road salting, wastewater) or when chloride is present in the bedrock (Ping et al., 2014). Then, preferential flows have been regularly highlighted as an error source since the associated water fluxes are not well sampled by zero-tension lysimeters (Tyler and Walker, 1994; Nkotagu, 1996). Finally, this method displays better results when rainfall and soil water are richer in chloride (e.g. sites close to the sea with high marine deposits or with low drainage flux) because the chemical analyses are more accurate for higher concentrations (Sammis et al., 1982; Grismer et al., 2000).

On the other hand, modelling errors could explain the bias presence. One of them could be the overestimation of the transpired water amount. However, deep drainage tends to produce during winter, while transpiration only takes place during the vegetation period (spring and summer). Therefore, if transpiration was overestimated we should observe an underestimation of the EW during spring and summer (low values), which is not the case (Fig. 7).

Hanson et al. (2004) measured deep drainage at the watershed level by accounting for rainfall and stream flow outputs and compared their measurements with the predictions of several models. Their multimodel comparison displayed similar RMSE (65.5 to 225.6 mm) and relative bias (−27.6 % to 20.5 %) values, but the Pearson's coefficient displayed by HETEROFOR is definitely located in the high tail of the study range (0.61 to 0.95). However, the performances of their models are not strictly comparable to ours since the reference method for estimating drainage differs (Sammis et al., 1982; Grismer et al., 2000; Obiefuna and Orazulike, 2011).

4.3 Simulating phenology and water balance in heterogeneous stands

Increasing the functional trait diversity and promoting uneven-aged stands are among the management strategies that foresters can use to make their forests more resistant to stressful conditions and more resilient after a disturbance (Pedro et al., 2015; Jactel et al., 2017; Anderegg et al., 2018). With the growing interest for mixed and uneven-aged stands, various attempts have been made to better account for stand structure in process-based forest models. Some of these models present very detailed 3D representations of individual tree structure but describe generally only specific physiological processes (e.g. LIGNUM, EMILION and MAESPA). Such models are very useful tools for analysing outcomes of ecophysiological experiments and obtaining a better understanding of specific ecophysiological processes (e.g. drought sensitivity) in structurally complex stands. Since they are generally computationally expensive and applied to one or a limited number of individuals, they can however not be used for simulating long-term forest dynamics according to various climate and forest management scenarios. Other individual-based models can be applied on all the trees of a stand in long-term simulations but at the cost of a coarse representation of physiological processes (e.g. SORTIE/BC). These models are interesting for the analysis of tree growth dynamics in heterogeneous forests but are less suitable for taking into account the changing environment. Since they simplify stand structure representation, cohort-based models can afford a detailed process-based description of the main processes involved in tree growth (e.g. 4C, ANAFORE, PSIM-DNDC and 3D-CMCC; see Table 5 for model characteristics). Here, the compromise is made on the spatial representation, which accounts for the vertical gradient in growing conditions but not for the horizontal heterogeneity. Such models can indeed not distinguish between stands composed of the same trees but with various degrees of spatial aggregation (e.g. intimate vs. patch-wise mixture). Similarly, some individual-based models choose to sacrifice the horizontal heterogeneity of some processes (e.g. iLand and Hybrid, which calculate most of the water balance at stand scale; see Table 5 for model characteristics).

To simulate the impact of management in heterogeneous forests under changing conditions, we developed a spatially explicit individual-based approach designed to account for individual variability, local environment and adaptive behaviour of trees (Berger et al., 2008). The compromise was not achieved by strongly reducing the complexity of a particular aspect (spatial representation, process description or spatial or temporal coverage), but instead we tried to develop a balanced approach in which each aspect is described with the same level of complexity.

Among the existing individual-based models, BALANCE and NOTG-3D are close to HETEROFOR since they were designed according to the same philosophy. They present however some substantial conceptual differences (Table 5). Except BALANCE for leaf yellowing, HETEROFOR is the only model determining budburst, leaf yellowing and fall at the tree level. While rainfall partitioning is only calculated in HETEROFOR, the spatial representation of local climate conditions in the canopy is finer in BALANCE and NOTG-3D, which consider different canopy layers or voxels. Regarding transpiration, HETEROFOR and BALANCE implement the widely used Penman–Monteith equation, while it is determined as part of detailed energy budget in NOTG-3D. Finally, they all describe soil water dynamics at the individual scale, but HETEROFOR displays a more mechanistic approach for describing soil water transfer among horizons (bucket vs. Darcy model).

In HETEROFOR, some processes were described at two spatial scales (tree or stand level) in order to have the opportunity to compare the two approaches and choose the most appropriate one depending on the pursued objective. The phenological timing is species dependent and can optionally vary with tree size. This option (phenology at tree level) is very interesting since it accounts for both the ontogeny effect and the vertical gradient in climate conditions. With this option, a longer vegetation period is assigned to the understory compared to the overstorey, which allows for the improvement of radial growth predictions by correcting the growth underestimation in small trees and the overestimation in bigger ones (Appendix D). This first attempt to describe phenology at tree scale is quite empirical and could be adapted in the future as knowledge on inter-individual phenology differences improves. Individual phenology observations for trees of all social status will be necessary to better calibrate and evaluate this module in an iterative cycle of model improvement.

For the water balance, HETEROFOR accounts for a direct tree size effect on stomatal conductance (stomatal conductance is inversely proportional to the height of largest crow extension) and for an indirect effect on the sapwood to leaf area ratio whose components both depend on tree size (Jonard et al., 2020a). In addition, individual transpiration is a function of the radiation intercepted by the tree, the local wind speed and the soil water availability. Finally, the tree adaptive behaviour to the local environment is described by an adaptation of the foliage biomass to local competition conditions and by specific leaf area varying with crown position within the canopy (Jonard et al., 2020a). Whatever the considered scale (tree or stand), HETEROFOR was able to correctly reproduce individual tree transpiration. Additional sap flux measurements as well as a characterization of the horizontal soil water content heterogeneity (using ground-penetrating radar tomography – GPR – technique for example) would be very useful to further evaluate the model performances and still enhance its ability to describe the complex hydrological functioning of heterogeneous forest. Among the possible improvements, mortality representation could be enhanced by considering hydraulic failure during severe droughts (Martin-StPaul et al., 2017). Another model improvement would be to take the interaction between the water cycle and the phenology into account by integrating a drought effect on budburst, leaf yellowing and fall as reported in some observation studies (Sanz-Pérez and Castro-Díez, 2010; Xie et al., 2018).

5 Conclusion

In this paper, two key modules of HETEROFOR are described in detail and evaluated in four sites/six stands. The phenological module correctly predicts the leaved period, which is essential to simulating light interception by trees, evapotranspiration, photosynthesis and respiration. With the hydrological module, HETEROFOR properly estimates rainfall interception, individual transpiration, soil water and deep drainage. Reproducing correctly the soil water dynamics is necessary to adequately predict photosynthesis since stomatal conductance closely depends on it. In addition, the description of the nutrient cycling requires accurate estimates of the water fluxes since water is the main vehicle for nutrient transport.

Our spatially explicit individual-based approach allows for a description of the phenology and water balance in structurally complex stands by partly accounting for the tree size effect on phenology and on tree transpiration, for the local environment modification (radiation and water availability) and for the adaptive behaviour of trees to local conditions (e.g. tree leaf area). Given the complexity of the functioning of heterogeneous forests, there are still a lot of ways to explore for improving the model, which will be done progressively as part of an iterative approach based on the comparison of predictions with observations. Our model will also be used to compare various spatial representation scales (tree, cohort, stand) and to determine the most appropriate one depending on the considered process and the pursued objective.

Simulating resource availability properly is necessary to produce robust predictions of tree growth under changing climate conditions. The next steps will be to extend the model validation to other European sites in order to cover a larger range of ecological conditions and to use HETEROFOR to simulate stands dynamics under various management options and climate scenarios.

Appendix A

Table A1 Description of some model parameters for European hornbeam (regarding light interception, respiration, carbon allocation and tree dimension increment) and origin of their value.

Appendix B

Figure B1 Comparison of the log-transformed observed and predicted monthly throughfall in Chimay (sessile oak), Louvain-la-Neuve (European beech) and Virton (European beech) between 2000 and 2016. The shaded area represents the confidence interval of the Deming regression (95 %) of observations on predictions, and the solid line corresponds to the identity line.

Appendix C

Figure C1 Temporal dynamics of soil-extractable water simulated with the tree approach in the three stands of Baileux for 2003. The shaded area represents the 80 % confidence interval of the values obtained for the various pedons. For comparison, the mean extractable water calculated with the stand approach is represented with a dashed line.

Appendix D

Figure D1 Comparison of observed and predicted basal area increments for sessile oak and European beech considering the two phenology modalities (tree vs. stand scale). The quality of predictions is indicated by the Pearson's correlation coefficient (r). The shaded area represents the confidence interval of the Deming regression (95 %) of observations on predictions, and the solid line corresponds to the identity line.