Radiative transfer modeling in structurally complex stands: towards a better understanding of parametrization

Abstract

• Key message

The best options to parametrize a radiative transfer model change according to the response variable used for fitting. To predict transmitted radiation, the turbid medium approach performs much better than the porous envelop, especially when accounting for the intra-specific variations in leaf area density but crown shape has limited effects. When fitting with tree growth data, the porous envelop approach combined with the more complex crown shape provides better results. When using a joint optimization with both variables, the better options are the turbid medium and the more detailed approach for describing crown shape and leaf area density.

• Context

Solar radiation transfer is a key process of tree growth dynamics in forest.

• Aims

Determining the best options to parametrize a forest radiative transfer model in heterogeneous oak and beech stands from Belgium.

• Methods

Calibration and evaluation of a forest radiative transfer module coupled to a spatially explicit tree growth model were repeated for different configuration options (i.e., turbid medium vs porous envelope to calculate light interception by trees, crown shapes of contrasting complexity to account for their asymmetry) and response variables used for fitting (transmitted radiation and/or tree growth data).

• Results

The turbid medium outperformed the porous envelope approach. The more complex crown shapes enabling to account for crown asymmetry improved performances when including growth data in the calibration.

• Conclusion

Our results provide insights on the options to select when parametrizing a forest radiative 3D-crown transfer model depending on the research or application objectives.

1 Introduction

Solar radiation provides energy to forest ecosystems and governs many biological and ecological processes such as photosynthesis, evapotranspiration, respiration and phenology. Its sharing among trees in the over- and understorey will influence tree regeneration, growth and survival, biomass allocation and tree morphology, as well as ground vegetation dynamics and diversity. Detailed knowledge and understanding of the radiation transfer within and beneath forest canopies is therefore of prime importance to quantify competition for light among trees and develop scientifically based sylviculture. Moreover, continuous-cover forestry and tree species diversification necessitate to maintain understorey light conditions favorable to the growth and survival of seedlings as well as to the mixing of species with contrasted shade tolerances.

However, accurate characterization of the light environment within forests through direct measurements is labor intensive, time consuming and nearly unrealistic routinely given the specific equipment and the large number of sensors required to capture the large spatial and temporal variability of the transmitted radiation (Baldocchi and Collineau 1994; Lieffers et al. 1999). Besides, indirect assessments of light availability are possible through canopy cover or canopy closure estimations by means of sensors or instruments (e.g., LAI 2000, spherical densiometer) or via hemispherical photographs (Jennings et al. 1999; Baudry et al. 2014; Fournier and Hall 2017; Hossain and Comeau 2019). Yet, these indirect methods also present shortcomings (Lieffers et al. 1999; Weiss et al. 2004; Chianucci and Cutini 2012; Yan et al. 2019). A workaround to these limitations is the use of models describing the radiative transfer through the forest canopy (Brunner 1998; Lieffers et al. 1999; Ligot et al. 2014a).

As reviewed by Ligot et al. (2014a) and summarized here, radiative transfer models may be classified according to the approach used to represent the canopy:

• (i) Models considering the canopy as one or several homogeneous horizontal layer(s), which makes these one-dimensional (1D) models particularly adapted to pure even-aged stands (Govind et al. 2013) though they were also applied to mixed stands (Kim et al. 2011);

• (ii) Three-dimensional (3D) crown models in which individual tree crowns, and sometimes trunks, are represented by spatialized geometric shapes (e.g., spheres, ellipsoids, cylinders, cones) and light interception by trees calculated with a ray tracing approach. The crown shapes can be either simple or composite through the assembly of several crown fractions for a more realistic representation, accounting notably for crown asymmetry or for foliage density variations within the crowns (Wang and Jarvis 1990; Gersonde et al. 2004). Compared to 1D models, they allow to consider stand heterogeneity and to predict the spatial variability of irradiance. However, their more detailed canopy representation necessitates a larger set of input data and parameters than 1D models;

• (iii) “Summary” light models with canopy representation intermediate to 1D and 3D models such as that proposed by Forrester (2014), allowing to account for canopy gaps and for multiple canopy components (i.e., age or dominance classes, species) through subdivision of the canopy into horizontal layers which may themselves consist of several components whose crowns overlap vertically;

• (iv) 3D-surface models using very detailed canopy mock-ups, representing individual tree organs with surfaces. Given the initialization data and parameter requirement, such models are generally applied to single trees or sparse canopies (Da Silva et al. 2008; Leroy et al. 2009).

Aside the canopy representation, models also differ depending on the approach used to determine the radiation attenuation when penetrating through the canopy:

• (i) The turbid medium formulation which follows the Beer-Lambert law and basically considers the canopy or its items (e.g., layer, crown, crown part) as consisting of small randomly scattered and oriented elements which do not reflect or transmit light.

• (ii) The porous envelope method in which canopy items are characterized by an “openness” empirical parameter specifying the probability of a ray to be intercepted.

These two approaches will be further described later (see Section 2.1.3).

When implementing, selecting or parametrizing a radiative transfer model, the modeler or the model user faces a series of questions regarding the choice of the model features. The first step concerns the canopy representation. Though selecting among the 1D, 3D-crown, summary light model and 3D-surface alternatives is rather easy based on forest type, on research or application objectives, and on available data, choosing from sub-alternatives such as, for instance, the kind of crown shape and in particular its complexity (i.e., simple vs composite) for 3D-crown models is less straightforward. The second step consists in choosing between the turbid medium and the porous envelop approaches to model the light interception by trees. Finally, the nature of the data used to calibrate the model and/or to evaluate its accuracy may also be crucial. Indeed, one may wonder if the same level of complexity is needed to predict the transmitted radiation reaching the ground than to estimate light interception by individual tree crowns.

This study aims to provide an efficient and robust framework to parametrize spatially explicit models of light interception. We selected for this purpose a 3D-crown light interception library (SAMSARALIGHT, http://hdl.handle.net/2268/187361) and used it coupled to HETEROFOR, a spatially explicit model describing individual tree growth based on resource sharing in heterogeneous forests (Jonard et al. 2020; de Wergifosse et al. 2020). To meet the objective of this study, the coupled model was calibrated for several options: crown shapes of contrasting complexity to account for their asymmetry, the turbid medium vs the porous envelope approaches. Moreover, the added value of accounting for the intra-specific variations in leaf area density was also evaluated for the turbid medium approach. Furthermore, while radiative transfer model calibration is usually performed based on direct or indirect measurements of transmitted radiation, we also used tree growth data as a proxy of the intercepted radiation by individual crowns and combined the two response variables in both sequential and joint optimization schemes. Model performances were evaluated for all calibration runs and were compared so as to identify the best model options depending on the considered response variable(s). In all cases, calibration was carried out by resorting to Bayesian optimization.

2 Materials and methods

2.1 HETEROFOR model

2.1.1 General overview

HETEROFOR is a spatially explicit, tree individual-based and process-based model describing forest stand dynamics based on resource use (light, water and nutrients) and silvicultural operations. It is hosted in the Computer-Aided Projection of Strategies In Silviculture (Capsis) platform (Dufour-Kowalski et al. 2012), dedicated to the simulation and the modeling of tree growth and stand dynamics. The overall operation of HETEROFOR, the modeling of the carbon-related processes (i.e., photosynthesis, respiration, carbon allocation, tree dimensional growth) and the water balance and phenology modules have been described in Jonard et al. (2020) and de Wergifosse et al. (2020) together with first validations of the model on these aspects. The present paper focuses on the description and the calibration of the light interception module.

2.1.2 Canopy structure modeling

HETEROFOR being spatially explicit and individual-based, each tree is positioned within the investigated stand and its aboveground architecture is represented on the basis of dendrometric variables either provided by inventories for the initialization step or derived from growth modeling for the simulation steps. Each trunk consists of a cylinder with circumference corresponding to the girth at breast height (C130, cm), and extending from the forest floor to the height of crown base (Hcb, m). Crowns may be depicted considering either single ellipsoid (shape “E”), two half-ellipsoids representing the lower and the upper parts of the crown (shape “B”), paraboloid (shape “P”) or conical (shape “C”) geometric representations. For each of these shapes, the crown may be either trunk-centered (“c”) or decentered (“d”). Furthermore, in order to account for the crown plasticity in response to competition, more complex crown forms considering eighths of ellipsoids (shape “M”), fourths of paraboloids (shape “Pm”) or fourths of cones (shape “Cm”) are also implemented, resulting in a total of 11 possible crown types in HETEROFOR (Fig. 1). In each case, the crown horizontal positioning and extension are determined from the crown radii in each of the four cardinal directions. For the trunk-centered crown types (i.e., “Ec,” “Bc,” “Pc” and “Cc”), the lengths of the south-to-north and of the west-to-east horizontal axes are defined as the sums of the crown radii in the corresponding directions and these axes are centered on the x- and y-trunk coordinates. The decentered crown types are positioned so that the ellipsoid(s) representing the “Ed” and “Bd” crowns or the horizontal cross-section ellipse representing the base of the “Pd” and “Cd” crowns pass through the four points defined by the crown radii. For the “M,” “Pm” and “Cm” forms, each crown fraction is generated from, respectively, an ellipsoid, a paraboloid or a cone centered on the trunk and with horizontal semi-axes corresponding to the two crown radii delimiting that fraction. Finally, the lengths of the vertical semi-axes of the ellipsoid fractions used to represent crowns with lower and upper parts (i.e., “Bc,” “Bd” and “M” types) are determined as the difference between the height of the largest crown lateral extension (Hlce, m) and Hcb for the lower part and as the difference between the total tree height (Htot, m) and Hlce for the upper part. For all other crown types, the crown length is set to the difference between Htot and Hcb. Given their development type (sympodial vs monopodial), crowns of broadleaved species are usually depicted with ellipsoidal crown shapes while paraboloid and conical shapes are mostly used to represent coniferous crowns.

Fig. 1

Geometrical representations of possible crown shapes in HETEROFOR

2.1.3 Radiative transfer modeling

The SAMSARALIGHT library used to model the transfer of solar radiation through the canopy is based on a radiative transfer model developed and implemented in the Capsis platform by Courbaud et al. (2003) which was later progressively improved (Ligot et al. 2014b). SAMSARALIGHT is particularly suitable for uneven-aged and mixed forests as it provides both the energy intercepted by each tree or tree part (i.e., crown, crown part and/or trunk) and the distribution of irradiance on the ground, while requesting a quite limited number of input parameters. In SAMSARALIGHT, the stand area is subdivided into cells arranged as a regular grid on the ground and the model describes how the energy of light beams pointing to the center of each of these cells attenuates as they penetrate into the canopy.

In SAMSARALIGHT, the first step is the creation of beams from the incident radiation measured at a meteorological station representative of the climate of the study site. The measured global radiation is partitioned into direct and diffuse energies, through coefficients either fixed to arbitrary values or derived from measurements or from empirical equations accounting for atmospheric clearness through the global-to-extraterrestrial radiation ratio as proposed by Erbs et al. (1982) or Jacovides et al. (2010). Direct and diffuse energies are then distributed among beams. Direct beams are generated for regular hour angle intervals along the mean solar path for each month of the growing season according to astronomic laws providing the sun height angle and the sun azimuth for each hour angle (Bonhomme 1993). Direct energy is then shared between the beams in proportion to the sinus of their height angle. Diffuse beams are created at regular azimuth and elevation intervals on the sky hemisphere and diffuse energy is shared between the beams according to either the “Uniform Overcast Sky” (UOC) or the “Standard Overcast Sky” (SOC) distributions (Moon and Spencer 1942). The former distribution assumes identical radiation flux from all sky directions while the latter considers variation as a function of height angle, with larger radiation flux from sky directions towards the zenith.

In a second step, each generated beam is pointed to the center of each cell and intersections with intercepted trees are calculated. Based on crown geometric representations as presented in Section 2.1.2, an entry point and an exit point are determined for each beam intersecting a crown (see Courbaud et al. (2003) and Ligot et al. (2014b) for computation details). Not considered in the initial version of the model, trunks may now be optionally represented as cylinders which, in contrast to crowns, do not transmit light and interrupt the path of beams.

In a third step, the energy intercepted by each tree and transmitted to the ground is calculated using either the turbid medium or the porous envelop approaches. The turbid medium approach assumes that radiation extinction inside a crown (or a crown part) follows that of a monochromatic ray through an homogeneous medium, which is described by the Beer-Lambert law. According to this law, light transmission is a function of the medium absorption coefficient and of the light path length inside the medium (l, m). In the case of a canopy, the absorption coefficient may be expressed as the product of an extinction coefficient (k, -) depending on the leaf and branch orientation and spatial distribution, leaf area density (LAD, m²m^− 3) and a clumping index (Ω, -) accounting for the aggregation of leaves and branches, leading to the following formulation:

$$ \frac{I_{out}}{I_{in}} = exp(-k \cdot {\Omega} \cdot LAD \cdot l) $$

(1)

where I_in and I_out are the irradiances at the entry point and at the exit point, respectively. The path length within the crown (or the crown part) is given by the distance between the beam entry and exit points. Two options are possible for determining LAD, either fixing it to a specific value (LAD_constant option) or modeling it by considering the specific leaf area (SLA, m² kg^− 1) variation depending on the vertical position within the canopy (LAD_model option). Based on observations and literature references (e.g., Fellner et al. 2016; Forrester et al. 2017; Leuschner and Meier 2018), SLA is assumed to decrease from the bottom of the canopy to its top. Therefore, when considering the LAD_model option, the crown SLA for a given tree is evaluated by linear interpolation of reference SLA values for the top (SLA_top, m² kg^− 1) and for the bottom (SLA_bottom, m² kg^− 1) of the canopy as a function of crown height level. The top and bottom of the canopy are determined based on the mean height of the lower/upper crown part of the n lowest and highest trees in the vicinity of the target tree, respectively. LAD is obtained by multiplying the SLA by the leaf biomass and by dividing it all by the crown volume derived by applying the ellipsoid, the paraboloid or the cone volume formula depending on the crown form. For crowns subdivided into lower and upper parts, the sharing of leaf biomass between both parts was taken into account through the upper fraction leaf biomass (UFLB) parameter (see Appendix for more details). The determination of the other parameters will be discussed below. In contrast to the turbid medium analogy, the porous envelop approach was formulated considering a fixed crown openness (p, -), not depending on the light path distance through the crown neither on foliage density or arrangement. In this case, radiation interception inside the crown (or the crown part) may be expressed as:

$$ \frac{I_{out}}{I_{in}} = p $$

(2)

When intercepted by a crown (or a crown part), the energy of a beam is decremented by the quantity determined according to the chosen approach and, in parallel, the energy intercepted by the corresponding crown (or crown part) is incremented by this quantity. This beam radiation extinction occurs for the successive trees traversed by the beam until it reaches the target cell. In case a beam hits a trunk, the trunk intercepted energy is incremented by the remaining energy of the beam and the beam energy is set to zero. After processing of all the beams for each cell, the total radiation reaching a ground cell is given by the sum of the remaining energies of the beams coming from all sky directions, and the total energy intercepted by a tree part (i.e., crown, crown part or trunk) equals the sum of the energies intercepted from the different beams. Besides this actual energy interception by trees, potential energy interception is also determined for each tree and corresponds to the energy that would have been intercepted by this tree in the absence of its neighbors. The ratio between actual and potential interception may be considered as a light competition index (LCI, -). For the sake of computing efficiency, the algorithm has been optimized in order to, for each beam, limit beam-tree interception computations only to trees likely to effectively intercept that beam. We refer to Courbaud et al. (2003) for details about this optimization technique. Furthermore, virtual wrapping of the plot as a torus is implemented in the model to avoid erroneous radiative balance around plot borders. It allows repetition of the plot in every direction and, thereby, avoids overestimation of incident radiation coming for the edges and provides exact radiation balance at the stand level (see Courbaud et al. (2003)). Finally, speedup of the algorithm was recently achieved through parallelization of the radiative balance computations for the different cells.

2.1.4 Photosynthesis and primary production

In HETEROFOR, annual gross primary production (GPP, kg of C year^− 1) may be simulated for each tree using either an empirical method based on photosynthetic active radiation (PAR) use efficiency or a process-based approach for modeling photosynthesis. In the present study, we adopted the process-based approach, which resorts to the CASTANEA library of Capsis (Dufrêne et al. 2005) in which the biochemical model of Farquhar et al. (1980) is implemented. The Farquhar’s model notably requires, for each tree, the direct and the diffuse components of the photosynthetically active radiation (PAR) absorbed at an hourly time step per leaf area unit. These absorbed direct and diffuse PAR components are determined by multiplying the direct and diffuse energies intercepted by tree crowns provided by the SAMSARALIGHT library (see Section 2.1.3) by (i) the canopy reflection coefficient (set to 0.11) to account for albedo and (ii) the proportion of PAR in solar radiation (set to 0.46). The obtained quantities are then divided by the crown sunlit leaf area and by the crown total leaf area for, respectively, the direct and diffuse absorbed PAR. Furthermore, stomatal conductance of the canopy is modulated both by air relative humidity and by soil water potential based on the approach of Ball et al. (1987) to account for the stomatal response to these variables (see Jonard et al. (2020) for details). Simulated GPP is then converted to individual annual net primary production (NPP, kg of C year^− 1) by considering a ratio between NPP and GPP (r_NppGpp, -) which accounts for carbon losses due to maintenance and growth respiration.

The NPP to GPP ratio is acknowledged to vary with tree characteristics (Mäkelä and Valentine 2001) and preliminary analyses based on our tree growth data (see Section 2.2.2) revealed that its variability was generally largely explained by the crown to stem diameter ratio (Dd) according to a linear relationship. The Dd ratio was standardized in order to remove the effect of tree size on its variation and to provide thereby an index reflecting tree past competition conditions, leading to the following equation:

$$ r_{NppGpp} = \alpha_{r_{NppGpp}} + \beta_{r_{NppGpp}} \cdot Dd_{Index} $$

(3)

where $\alpha _{r_{NppGpp}}$ and $\beta _{r_{NppGpp}}$ are parameters and Dd_Index (-) is defined as

$$ Dd_{Index} = \frac{Dd}{Dd_{pred}} $$

(4)

with Dd_pred (-) being the crown to stem diameter ratio predicted based on C130:

$$ Dd_{pred} = \alpha_{Dd_{Index}} + \beta_{Dd_{Index}} \cdot C130 + \gamma_{Dd_{Index}} \cdot \frac{1}{C130} + \delta_{Dd_{Index}} \cdot \frac{1}{C130^{2}} $$

(5)

in which $\alpha _{Dd_{Index}}$, $\beta _{Dd_{Index}}$, $\gamma _{Dd_{Index}}$ and $\delta _{Dd_{Index}}$ are parameters.

As described by Jonard et al. (2020), the carbon allocated each year to the growth of structural tree compartments (Δb_structural, kg of C year^− 1) is then determined as follows:

$$ {\Delta} b_{structural} = NPP + rt_{leaf} + rt_{fine \ root} - p_{leaf} - p_{fine \ root} - p_{fruit} - s_{branch} - s_{root} $$

(6)

where rt_leaf and rt_{fine root} are the annual leaf and fine root retranslocations, p_leaf,p_{fine root} and p_fruit are the annual leaf, fine root and fruit productions, and s_branch and s_root are the carbon fluxes compensating for annual branch and root mortality. All terms of this equation are expressed in kg of C year^− 1 and are evaluated using coefficients or parameters either determined from literature or fitted on field observations.

Subsequently, the increment in structural biomass (Δb_structural) is partitioned between above- and below-ground tree parts, Δb_{structural above} and Δb_{structural below} (kg of C year^− 1), respectively. The combined simulated increment in DBH and Htot (ΔD²H_sim, m³ year^− 1) is then derived from Δb_{structural above} using:

$$ {\Delta} D^{2}H_{sim} = \frac{\Delta b_{structural \ above}}{\beta_{woody \ above} \cdot \gamma_{woody \ above} (D^{2}H)^{\gamma_{woody \ above}-1}} $$

(7)

where β_{woody above} and γ_{woody above} are parameters from the following allometric equation (e.g., Genet et al. 2011; Hounzandji et al. 2015 ):

$$ b_{\textit{structural above}} = \alpha_{\textit{woody above}} \cdot \ \beta_{\textit{woody above}}(D^{2}H)^{\gamma_{\textit{woody above}}} $$

(8)

Simulated tree height increment (ΔHtot, m year^− 1) and trunk diameter increment (ΔDBH, cm year^− 1) may then be determined by developing and rearranging ΔD²H_sim, and simulated increments in root, stem and branch biomasses are obtained following allocation rules based on allometric relationships. We refer to Jonard et al. (2020) for a complete description of primary production and growth modeling.

In other respects, ΔD²H may be easily determined from repeated stand inventory data, providing thereby observed values for this variable (see Section 2.2.2). These ΔD²H_obs values will be compared to ΔD²H_sim in order to estimate model parameters through optimization (see Section 2.3.1).

2.2 Measurements

HETEROFOR calibration was carried out by resorting to two data sets, namely, hemispherical photographs and repeated stand inventories. Hemispherical photographs allowed assessing undercover radiation and its spatiotemporal variation while repeated stand inventories permitted to quantify individual tree growth. These quantities were then compared to model outputs so as to estimate parameters of the radiative transfer and of the growth modules. The two next sections present each data set in details. The methods used to estimate parameters and to evaluate model performances will be described subsequently (see Section 2.3).

2.2.1 Understorey radiation with hemispherical photographs

We benefited from hemispherical photographs performed in July 2010 by Ligot et al. (2013) within regeneration plots of 19 uneven-aged forest stands in the Belgian Ardennes dominated by European beech (Fagus sylvatica L.) and sessile oak (Quercus petraea (Matt.) Liebl.), covering a wide range of proportions of both species from pure beech stands to stands mainly dominated by oak. The photographs were acquired before sunrise with a Nikon D90 camera equipped with a Sigma 4.5 mm fisheye installed above the regeneration using a self-leveling mount. Between 6 and 39 photographs were taken in each stand, depending on the experimental site area. The photographs were processed using the PiafPhotem (Adam et al. 2006) and the GLA (Frazer et al. 1999) softwares so as to determine three indexes of light availability for the growing season (assumed as extending from 1st April to 31st October): the percentages of direct above canopy light (PACL_direct, %), of diffuse above canopy light (PACL_diffuse, %) and of total above canopy light (PACL_total, %). For validation of the technique, PACL_total estimates were successfully compared with PAR measurements carried out with sensors in five sites during one day in July 2010. Besides, in order to allow for simulations with SAMSARALIGHT, the collection of the hemispherical photographs was accompanied by dendrometric measurements and mapping of trees on concentric circular plots centered on each camera location: all trees with C130 larger than 20 cm were considered on a 15-m-radius plot and all trees with C130 larger than 7.5 cm were inventoried on a 7-m-radius plot. We refer to Ligot et al. (2013) and Ligot et al. (2014b) for a more complete description of the investigated sites as well as for more details on the processing of the hemispherical photographs and its validation. This 2010 data set was completed with new hemispherical photographs and dendrometric measurements carried out by the same authors on 4 additional regeneration experimental sites in 2011 as well as with repetitions of the measurements for subsets of the 19 initial stands in 2011 and 2012. Furthermore, hemispherical photographs were also collected in 2012 within the three plots of the Long-Term Ecosystem Research (LTER) site of Baileux (see Section 2.2.2). This resulted in a total of 654 hemispherical photographs (i.e., 1 962 individual PACL_direct, PACL_diffuse and PACL_total values) characterized by their local canopy environment. These data will be referred to as the “radiation data set” in the following, they are available through an open-access repository (André et al. 2021).

A large variability of understorey light conditions was encountered among the measurement locations, with percentage of above canopy light from 1% to 61% in oak-dominated plots and from 2% to 43% in beech-dominated plots, reflecting the wide gradient of canopy cover ranging from a close canopy to a large gap size of about 1 200 m² without overstorey.

The radiative balance was run on each site and for each measurement year, providing simulated PACL_direct, PACL_diffuse and PACL_total values in order to assess model parameters through optimization (see Section 2.3).

2.2.2 Inventories and individual tree growth

Stand inventory data from three Belgian level II plots (i.e., Louvain-la-Neuve, Chimay and Virton sites) of the International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests (ICP-Forests) and from the LTER site of Baileux were used to provide growth measurements at the individual level. The ICP-Forests level II plots consist each of a 0.25 ha (50 × 50 m) central zone surrounded by a 15 m wide buffer zone. The Louvain-la-Neuve plot is installed in an even-aged nearly pure beech stand growing on a leached brown soil (Abruptic Luvisol), the site of Chimay is located in a sessile oak and hornbeam (Carpinus betulus L.) coppice-with-standards stand on an acid brown soil characterized by the presence of a clay layer at a depth ranging between 40 and 80 cm (Dystric Cambisol) while the Virton plot consists of a mixed hardwood stand dominated by beech with ash (Fraxinus excelsior L), pedunculate oak (Quercus robur L.) and sycamore maple (Acer pseudoplatanus L.) as the main secondary species on a well-drained brown soil (Calcaric Cambisol). The Baileux LTER site is located in a 60 ha mixed sessile oak and beech forest lying on acid brown earth soil (Cambisol) on a plateau in the western part of the Belgian Ardennes. Three plots with central and buffer zones were installed at this site in stands dominated either by oak or by beech as well as in a balanced mixture of both species. Beyond oak and some beech trees, the oak-dominated plot also presents a hornbeam coppice in the understorey. We refer to Jonard et al. (2020) and to de Wergifosse et al. (2020) for a more complete description of these study sites.

In each of these plots, all the trees with C130 larger than or equals to 15 cm were spatialized and inventoried with the identification of the species and the measurement of C130, Htot, Hcb, Hlce and crown radii. A first inventory was carried out during winter 1999–2000 for each ICP-Forests level II plot while it was performed during winter 2001–2002 for the Baileux site. Inventories were repeated by measuring the same variables 10 years later at each site, namely, during winter 2009–2010 for the ICP-Forests sites and during winter 2011–2012 for Baileux. In addition, litterfall data collected yearly on the ICP-Forests plots allowed quantifying tree fruit production specifically for each year on each of these sites. Such data were not available for Baileux and fruit production at this site was estimated for each species using an average allometric relationship fitted on all fruit litterfall data from the three studied ICP-Forests plots (see Jonard et al. (2020)).

Inventory data comparison allowed to determine the average growth observed for each tree over the inter-inventory period by expressing it under the variable ΔD²H_obs (m³ year^− 1) integrating increments for both DBH and Htot, defined as:

$$ {\Delta} D^{2}H_{obs} = \frac{DBH_{final}^{2}Htot_{final} - DBH_{initial}^{2}Htot_{initial}}{N_{years}} $$

(9)

where subscripts “initial” and “final” refer, respectively, to the initial and final inventories and N_years is the number of years between them.

For each study plot, the HETEROFOR model was applied considering the initial inventory as starting stand stage, which provided ΔD²H_sim values (see Eq. (7)). These growth predictions were compared to the observations determined by Eq. (9) so as to estimate model parameters. For computing efficiency in the optimization process, tree growth was simulated on a single year for which meteorological conditions were found representative of the average annual meteorological pattern encountered on each site over the inter-inventory period. The year 2004 was selected for this purpose. For the study sites, mean temperature and total rainfall depth over that year ranged between 9.6 and 11.2 ^∘C and between 885 and 1 156 mm, respectively.

As for the radiation data set, these growth data are freely available (André et al. 2021).

2.3 Model calibration and evaluation

2.3.1 Parameter optimization

Model calibration was carried out in order to provide estimates for the parameters involved in the radiative balance (Eqs. (1)–(2), and (A.2)–(A.11) in the Appendix) and in the growth (Eqs. (3)–(5)) modules. This calibration exercise focused on the three main species found at the study sites, namely, sessile oak, European beech and common hornbeam, and parameters were estimated separately for each of the three species. So as to simplify the optimization process and to avoid instability and non-uniqueness issues, the values for some of the parameters were determined either based on literature or through separate fitting.

Regarding the “Turbid medium” mode of the radiative balance, and more specifically the “LAD_model” option, parameters α_leaf, β_leaf and γ_leaf from Eq. (A.8) were set to values provided by Jonard et al. (2006) for oak and beech and were fitted using unpublished data for hornbeam, while values for parameters SLA_bottom and SLA_top were evaluated from Forrester et al. (2017). For the “LAD_constant” option, LAD was fixed for each species to the corresponding average value obtained from the “LAD_model” option. The crown (or crown part) volume (V ) and the light path length (l) are geometrical parameters, directly determined by the crown shape and dimensions and, for l, by the beam path. Among the remaining parameters, k and Ω are both factors of a same product and cannot be estimated separately without fixing one of them. As a result, Ω was set to 1, considering then no clumping of canopy elements within tree crowns. Therefore, k and UFLB, for the “Bc,” “Bd” and “M” crown shapes (“LAD_model” option), are the two optimized parameters for the “Turbid medium” mode.

In the case of the “Porous envelop” mode, the canopy openness (p) is the only free parameter and has been optimized.

Finally, parameters $\alpha _{r_{NppGpp}}$ and $\beta _{r_{NppGpp}}$ were both optimized to provide estimates of r_NppGpp, while Eq. (5) was fitted beforehand on data from the same study sites.

The parameters were optimized resorting to a Bayesian procedure with iterative Markov Chain Monte Carlo (MCMC) carried out through the DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm, using more specifically the DREAM_(zs) version for greater efficiency (Vrugt 2016). For the present study, the R version of the algorithm implemented within the “BayesianTools” package (Hartig et al. 2019) was used.

Two different optimization strategies were tested and compared. The first strategy consists in a two steps optimization, optimizing first the radiative balance parameters based on the radiation data and fixing then these parameters to the median value of their retrieved posterior distribution for the optimization of the r_NppGpp ratio parameters based on the stand inventory comparison data. The first and the second steps of this optimization strategy will be referred to as, respectively, “Radiation data only” and “Sequential Radiation → Growth” optimization schemes in the following. The second strategy combined the optimization of both the radiative balance and the r_NppGpp ratio parameters using both sets of data in one single run, it will be denominated as “Joint Radiation × Growth” optimization scheme. These optimization schemes were repeated for each crown form representative of the shape of broadleaved species (i.e., “Ec,” “Ed,” “Bc,” “Bd” and “M”), the three species of interest in the present study belonging to this group. Three Markov chains of 20 000 iterations were run in each case, starting from random initial parameter values drawn within their acceptable ranges. The burn-in period (i.e., the number of initial iterations required for the chains to reach convergence) was delimited based on visual inspection of the chain history plots and was set to 6 000 iterations for “Radiation data only” runs and to 10 000 iterations for “Sequential Radiation → Growth” and “Joint Radiation × Growth” runs. These initial iterations were discarded. Convergence was further checked based on the Gelman-Rubin statistical indicators (Gelman and Rubin 1992).

2.3.2 Statistics

Comparison of model fits for the different considered model options and crown shapes was carried out based on Bayes factors (Kass and Raftery 1995). The Bayes factor (BF) can be computed as the ratio between the marginal likelihoods of two compared models and is an indicator of relative evidence of one model over the other given the data used to fit the models. The degree of evidence may be scaled based on the BF value. Notably, considering the largest marginal likelihood as the numerator of the ratio, substantial evidence in favor of the corresponding model over the model associated with the marginal likelihood in the denominator may be considered for BF≥ 3.2 (or log₁₀(BF) ≥ 0.5). This pairwise comparison method may be extended to the comparison of several models in order to identify the most evidenced model(s) in one step. For this purpose, the model presenting the maximum marginal likelihood value among all compared models was considered as the reference (M^∗) and the Bayes factors were computed relatively to it, considering this maximum marginal likelihood value as the numerator and the marginal likelihood value of the alternative model (M_i) as the denominator:

$$ BF_{i} = \frac{\textrm{pr(\textbf{D}}|M^{*})}{\textrm{pr(\textbf{D}}|M_{i})} $$

(10)

in which pr(D|M) corresponds to the marginal likelihood of the data D under the hypothesis of model M and is obtained by:

$$ \textrm{pr(\textbf{D}}|M) = \int \textrm{pr(\textbf{D}}|\theta, M) \pi(\theta, M) d\theta $$

(11)

where 𝜃 is the optimized parameter vector, π(𝜃,M) is the parameter prior density and pr(D|𝜃,M) is the likelihood function of 𝜃. We refer to Kass and Raftery (1995) for more details on Bayes factor theory and computation.

Thereby, using Bayes factors, the compared models could be ranked in terms of evidential strength with regard to the data. It is worth noting that Bayes factors only allow to compare models fitted on a same data set.

Besides, model fits were also evaluated by examining the correspondence between model predictions and observations based on statistical tests and measures of agreement. Statistical tests consisted in a paired Student’s t test and in a linear regression between observations and predictions. The selected measures of agreement were the Pearson’s correlation coefficient (r), the root-mean-square error (RMSE) and the fractional bias (FB) (Janssen and Heuberger 1995). The paired Student’s t test and FB quantify the bias between compared values. The regression analysis tests for the significance of the deviation of the intercept and of the slope from 0 and 1, respectively. Deming regression was used to account for the fact that both the dependent and the independent regression variables (i.e., observations and predictions) were associated with errors. Finally, the correlation coefficient quantifies the strength of the linear relationship between the two compared variables and the RMSE expresses the discrepancy between paired values.

Posterior distributions of optimized parameters were summarized through their median value and their highest density credible intervals determined using the “hdi” function of the “bayestestR” package (Makowski et al. 2019). Credible intervals were established at the 90% level rather than at the generally more commonly used 95% level for computational stability reason (Kruschke 2014).

3 Results

Preliminary analyses of the radiation data revealed that, for a large number of measurement points, the radiation transmitted to the ground was mainly from beams which had not been intercepted by any tree crown. Consequently, the corresponding observations did not bring relevant information for the present study which focuses on the estimation of the parameters governing radiation interception by trees. Therefore, the measurement points for which more than 40% of the transmitted radiation came from unobstructed beams were discarded. Only the remaining data set will be considered in the following (Fig. 2).

Fig. 2

Proportion of the total radiation reaching the understorey from unobstructed beams as a function of the observed percentage of total above canopy radiation

For the “Radiation data only” optimization scheme (Fig. 3a), the turbid medium approach combined with the LAD_model option generally presents the lowest BF values whatever the crown form. For this configuration, the “Ed” crown form is associated with the minimum BF value (i.e., log₁₀(BF) = 0, corresponding to the “best model fit,” with the highest marginal likelihood for this optimization scheme) followed by crown form “Bd.” The other three crown shapes show similar BF values (log₁₀(BF) ≈ 1.5), above the threshold value indicating substantial evidence in favor of the “best model fit” (represented by the dashed horizontal line on the Figure). In contrast, the highest BF values are systematically observed when considering the porous envelop approach, with particularly high values for the “Bc,” “Bd” and “M” crown forms. Finally, BF values for the model configuration combining the turbid medium with the LAD_constant option are in general intermediate to those of the two preceding configurations, though sometimes close to (crown form “M”) or lower (crown form “Bc”) than that of the model configured with the LAD_model option. It is worth noting that the LAD_constant option often resulted in abnormal stand LAI values as LAD is not modulated by tree size and competition, thereby giving rise to unrealistic simulation outputs. Consequently, it was only presented here to provide an overview of the quality of its fits compared with the other investigated model configurations and it will not be considered in the following. Regarding the “Sequential Radiation → Growth” optimization scheme, BF values are systematically much lower when resorting to the porous envelop compared with the “turbid medium × LAD_model” model configuration (Fig. 3b). Both configurations show similar patterns for the evolution of BF values as a function of the crown form, with a general trend to increase from the simpler to the more detailed crown representations for types “E” and “B” while the “M” crown type presents the minimum BF value in each case. Therefore, the “M” crown form × porous envelop combination corresponds to the best model fit for this optimization scheme. On the other hand, for the “Joint Radiation × Growth” optimization scheme, lower BF values are generally found for the “turbid medium × LAD_model” model configuration compared with the porous envelop, exceptions being the “Ec” crown form for which the opposite is observed and the “Ed” crown form for which similar BF values are obtained for both model configurations (Fig. 3c). For this optimization scheme, the best model fit corresponds to the “M” crown form combined with “turbid medium × LAD_model” model configuration.

Fig. 3

Bayes factors of the models fitted for each crown form and each considered option according to (a) the “Radiation data only,” (b) the “Sequential Radiation → Growth” and (c) the “Joint Radiation × Growth” optimization schemes. For each “crown form × model configuration × optimization scheme” combination, the Bayes factor is determined as the ratio between the highest marginal likelihood value observed among all model fits carried out for the corresponding optimization scheme and the marginal likelihood value of the model fit for the given combination. The Bayes factors are presented in logarithmic scale for legibility of the figures. The dashed horizontal lines indicate the value 0.5, corresponding to the threshold above which substantial evidence appears in favor of the model fit with the highest marginal likelihood value (i.e., the model fit for which log₁₀(BF) = 0)

In the following, only the best model and optimization configurations will be further investigated, namely, the turbid medium with “Ed” and “Bd” crown forms for the “Radiation data only” optimization scheme (Fig. 3a), the porous envelop with the “M” crown form both for the “Sequential Radiation → Growth” scheme (Fig. 3b) and the turbid medium with the “M” crown form for the “Joint Radiation × Growth” scheme (Fig. 3c). In addition, the turbid medium with the “M” crown form for the sequential optimization will also be examined in more details as it provides comparable fit quality when considering together results for radiation (Fig. 3a) and growth (Fig. 3b). Statistical tests and measures of agreement proposed above to compare model predictions to observations are presented in Table 1 for these five top ranked configurations.

Table 1 Statistical comparisons between observations and model predictions

Focusing first on the comparison between predicted and observed transmitted radiations, one can notice that the three crown forms used with the turbid medium approach (“Ed,” “Bd” and “M”) in the “Radiation data only” or sequential optimization schemes show very similar results for all considered statistics. The diffuse radiation component always presents the FB values closest to zero and largely non-significant paired t test results while, for direct radiation, FB tends to indicate underestimation of model predictions and t tests are systematically significant. Intermediate results are found for total radiation. The lowest RMSE values, around 6.21%, are observed for total radiation and it increases slightly for diffuse and direct radiation with values around 6.67% and 6.99%, respectively. Values of correlation coefficient r are around 0.6 in all cases, though somewhat higher and lower for direct and diffuse radiation, respectively. Deming regressions reveal more contrasted behavior among radiation components, with intercept and slope values always around 0 and 1, respectively, for direct radiation while regression lines deviate from the 1:1 line for diffuse and, to a lesser extent, total radiations due to a tendency to underestimate observed values at intermediate radiation levels, and overestimated them at higher levels (Fig. 4). As already noticeable by examining Bayes factors (see Fig. 3a), much poorer correspondence between observations and predictions is observed with the porous envelop approach, for which higher and much lower values are respectively found for RMSE and r coefficients for all radiation components. This is corroborated by both the strong divergence of the Deming regressions from the 1:1 line and the large variability of their coefficients. Yet, the FB values slightly closer to zero and the larger P-values of the t tests observed for the porous envelop indicate a somewhat smaller general bias. The joint optimization strategy considering the “M” crown form with the turbid medium provides a quality of agreement among observations and predictions in-between that obtained for the turbid medium and the porous envelop with the “Radiation data only” or sequential optimization schemes. For this model configuration and optimization scheme, FB and t-ratio statistics indicate an overestimation of model predictions which is higher for diffuse than for direct radiation. Besides, this optimization configuration shows the largest RMSE values, while r values around 0.47 are intermediate. Finally, the slopes of the Deming regressions diverge more from the identity line compared with the turbid medium for the same crown form.

Fig. 4

Comparison of observed and predicted percentages of above canopy light (PACL) for (a, b, c) direct and (d, e, f) diffuse radiation considering the “Radiation data only” optimization strategy for the “Ed,” “Bd” and “M” crown forms, using the turbid medium radiative transfer formulation. The star symbols represent individual observed PACL values for each considered measurement location as a function of corresponding median values for the distributions of PACL predictions obtained through the application to the model of 1 000 parameter sample sets drawn from the posterior distribution. The dashed lines represent the Deming regressions on these couples of values and the dot-dashed curves delimit the 95% confidence intervals to the regressions. The error bars are the 90% prediction credible intervals. The solid lines correspond to the 1:1 relationships

With regard to the comparison between observed and simulated ΔD²H, the joint optimization with the turbid medium option and the “M” crown form shows the lowest t-ratios and the FB values closest to zero among the three selected configurations including optimization for growth. On the other hand, RMSE and r present similar values in all three cases, respectively around 0.02 m³ and 0.89. As illustrated in Fig. 5, Deming regressions are close to the 1:1 line when considering the turbid medium for both optimization strategies, while the slope tends to diverge towards large values for the sequential optimization with the porous envelop.

Fig. 5

Comparison of observed and predicted ΔD²H for (a) the “Joint Turbid medium × Growth,” (b) the “Sequential Turbid medium → Growth” and (c) the “Sequential Porous envelop → Growth” optimization strategies considering in each case the “M” crown form. The star symbols represent average observed ΔD²H values for each “species × C130 class × site” combination as a function of corresponding median values for the distributions of average predicted ΔD²H obtained through the application to the model of 1 000 parameter sample sets drawn from the posterior distribution. The dashed lines represent the Deming regressions on these couples of values and the dot-dashed curves delimit the 95% confidence intervals to the regressions. Vertical errors bars are the 90% confidence intervals to the observed ΔD²H average values, and horizontal error bars are 90% credible intervals for the predicted ΔD²H average values. The solid lines correspond to the 1:1 relationships

Parameter estimates of the radiative transfer model are presented in Fig. 6 for each of the three studied species. Median values of the extinction coefficient (k) around 0.32 and 1.26 are found for oak and hornbeam, respectively, without significant differences among crown forms nor among optimization schemes (Fig. 6a, c). Beech presents intermediate k values with, in this case, significantly lower values for the “Joint Radiation × Growth” than for the “Radiation data only” optimization schemes (Fig. 6b). Furthermore, for this latter species, k estimates are constantly around 0.47 whatever the crown form for the “Joint Radiation × Growth” optimization while they increase from around 0.75 for crown forms “Ec” and “Ed” to around 1.2 for crown forms “Bc,” “Bd” and “M” in the case of the “Radiation data only” scheme. Given the often large amplitude of their credible intervals, upper fraction leaf biomass (UFLB) estimates are generally not significantly different among species, crown forms and optimization schemes (Fig. 6d–f). The only exception is the “Radiation data only” optimization scheme for beech for which high UFLB values around 0.95 are observed for all crown forms while, excluding this case, median UFLB values amount to 0.46. Canopy openness (p) estimates of 0.34 are found for oak and beech while they amount to 0.47 for hornbeam (Fig. g–i). For each species, no significant difference is observed neither among crown forms nor optimization schemes, though somewhat lower p estimates are obtained for the “Radiation data only” compared to the joint optimization.

Fig. 6

Posterior estimates of (a, b, c) extinction coefficient (k), (d, e, f) upper fraction leaf biomass (UFLB) and (g, h, i) canopy openness (p) obtained for oak, beech and hornbeam for the different crown forms and for the “Radiation data only” and the “Joint Radiation × Growth” optimization schemes. Symbols represent median values of posterior estimate distributions and error bars correspond to parameter 90% credible intervals

Estimates for the NPP to GPP ratio are shown in Fig. 7 for each of the studied species, crown form, and optimization strategy (i.e., sequential or joint optimization) considering either the “turbid medium” or the “porous envelop” approach. In order to ease comparison of values among species and model options, estimates are directly presented for the r_NppGpp calculated based on the $\alpha _{r_{NppGpp}}$ and $\beta _{r_{NppGpp}}$ parameters with Eq. (3). So as to describe the amplitude of variation of r_NppGpp as a function of Dd_Index (Eq. (4)), estimates are provided for Dd_Index values spanning the index distribution, namely, the quantiles q_0.025, q_0.5 and q_0.975. For oak and hornbeam, r_NppGpp estimates decrease significantly as Dd_Index increases, ranging from values around 0.54 to 0.70 for oak and from around 0.35 to 0.45 for hornbeam. In contrast, for beech, no significant effect of Dd_Index on r_NppGpp is observed and estimate values around 0.47 are found. For oak, optimizing jointly the radiative balance and growth parameters considering the “porous envelop” approach provide, for all crown forms, systematically higher r_NppGpp estimates compared with the other optimization strategies. For low Dd_Index values of hornbeam, r_NppGpp estimates are systematically lower with the “turbid medium” approach than with the “porous envelop,” whatever the adopted optimization strategy (i.e., sequential or joint), though differences are not significant. In all other cases (i.e., low and intermediate Dd_Index for hornbeam, and beech), similar r_NppGpp values are retrieved whatever the optimization strategy. In addition, for a given species and optimization strategy, differences of r_NppGpp estimates among crown forms are never significant.

Fig. 7

Posterior estimates of NPP to GPP ratio retrieved for the different crown forms of each investigated species and for each optimization strategy. For oak and hornbeam, r_NppGpp estimates are presented for quantiles 2.5% (q_0.025), 50% (q_0.5) and 97.5% (q_0.975) of their respective Dd_Index distributions, while only average r_NppGpp estimates are shown for beech as no significant relationship between this parameter and Dd_Index was found for this species. Values for q_0.025(Dd_Index) are 0.7 and 1.0 for oak and hornbeam, respectively, q_0.5(Dd_Index) = 1.1 for the three species and q_0.975(Dd_Index) = 1.5 for both oak and hornbeam. Symbols represent median values of posterior estimate distributions and error bars correspond to 90% credible intervals

4 Discussion

Our results provide valuable information on the modeling approach and the options to select when setting up a forest radiative 3D-crown transfer model depending on the research or application objectives.

A first major outcome of the present study is the much lower quality of the fit found for the porous envelop approach compared with the turbid medium when performing the optimization on the radiation data set (see Fig. 3a), and the associated poor agreement between observed and predicted values (see Table 1). These observations contrast with literature reports indicating similar accuracy levels in model predictions for both approaches (Ligot et al. 2014a) or even sometimes better modeling performances for the porous envelop than for the turbid medium (Groot 2004). However, contrary to our parametrization assuming canopy openness as constant (see Eq. (2)), the latter author adapted this basic formulation which he denominated as a “hit model”. A first adaptation allowed to consider possible variation of canopy openness with the light exposure level, shaded crown parts being likely to present lower foliage density and, thereby, higher transmissivity than those exposed to high radiation levels. Secondly, canopy openness was expressed as a function of beam elevation angle, accounting in this way partly for the length of the beam path within the crown. These adapted formulations make the porous envelop closer to the turbid medium approach and explain the good performances reported by Groot (2004). Indeed, the turbid medium equation is explicitly path-length dependent (see Eq. (1)) and the LAD_model option allows for modulating the leaf density according to both tree dimensions and light competition conditions. Tree dimensions (C130, crown length and radii) are used to estimate leaf biomass (see Eq. (A.8) in Appendix) and crown (or crown part) volume in Eqs. (A.9)–(A.11), and SLA is dependent on the relative vertical position within the canopy. In this regard, Da Silva et al. (2012) show that the turbid medium and the porous envelop approaches are intrinsically linked when beam path length and foliage density are considered. In other respects, it is worth noting that the alternative proposed by Groot (2004) to the “hits model” was a “path-length model” which did not include LAD contrarily to the turbid medium formulation, thereby considering implicitly foliage density as constant and integrated within coefficient k. Therefore, it corresponds to our LAD_constant option which, beyond resulting in unrealistic simulated LAI values at some of our study stands, gave also rise to generally poorer fits compared with the LAD_model option (see Fig. 3a). Besides, Canham et al. (1994) argue that hit models would better fit than path-length models for species for which the foliage is mainly distributed at the periphery of the crown, and inversely for species showing more uniform distribution of leaves within the crown.

Focusing on the radiation data optimization using the turbid medium and the LAD_model option, the better fits obtained with the “Ed” and “Bd” crown types (see Fig. 3a) indicate the necessity to consider crown eccentricity for proper radiative transfer modeling within forest canopies. Nevertheless, the deterioration of the fit quality associated with the other crown forms for the same configuration appears as being rather minor, especially relatively to the much larger differences observed for the two other model formulations and options. These observations are consistent with findings from sensitivity analyses of forest radiative transfer models. According to these works, crown radius is among the most sensitive parameters for 3D crown models (Stadt and Lieffers 2000; Beaudet et al. 2002; Piboule 2001; Gersonde et al. 2004; Da Silva et al. 2012; Ligot et al. 2014a) while the crown shape plays a secondary role (e.g., West and Wells 1992; Larsen and Kershaw 1996; Brunner 1998; Piboule 2001; Piboule et al. 2005; Rojo et al. 2020). In our case, the crown representation is determined based on crown radius measurements for all modalities and slight differences appear only between the crown shape accounting or not for the crown asymmetry (“Ec,” “Bc” vs “Ed,” “Bd,” “M”). The crown shape effect reported in the literature is more marked than in our study since it was tested by comparing crown representations with contrasting geometrical features, while crown forms are all ellipsoids (or parts of ellipsoids) in our study.

In other respects, using the same Capsis library and partly the same radiation data set as in the present study (i.e., year 2010 only, without filtering for unobstructed beams), Ligot et al. (2014b) reported some overestimation of measured PACL_total for positions subject to high radiation levels as also observed in our results (see Table 1). These authors interpret these observations by the presence of small trees that might have affected PACL measurements based on the hemispherical photographs but which were not inventoried due to a stem circumference lower than the inventory threshold and, therefore, which were not accounted for in the model simulations. Moreover, the general poorer agreement found for PACL_diffuse compared with PACL_direct might also be explained by the presence of these small trees. Indeed, given their small size, these trees affect essentially light coming from the lower hemispherical angles, which represent permanently a part of the diffuse radiation while direct beams follow the course of the sun which is most of the time located at higher elevation angles. Yet, other authors also mentioned differences in fit quality and parameter sensitivity between direct and diffuse lights (Brunner 1998; Gersonde et al. 2004; Groot 2004). Finally, such overestimation of measured PACL might also partly arise from the fact that only leaves (i.e., no crown woody parts) are considered in the model for the computation of crown LAD (see Section 2.1.3) while shoots and branches are present on the hemispherical photographs used to derived the measurements. Furthermore, shoots and branches are likely to be more visible (i.e., less overlaid by foliage) on the photographs taken under light canopy cover, namely in zones more subject to high radiation levels such as those for which the largest PACL overestimations are found in our results.

The sequential optimization on growth data reveals a better fit quality for the porous envelop approach than for the turbid medium mode whatever the crown form (Fig. 3b). Nevertheless, these results for the porous envelop are associated with the poorest fits found for radiation data optimization (see Fig. 3a, Table 1). Moreover, the quality of the comparisons between observed and simulated ΔD²H values for the selected “top ranked” configurations does not allow to clearly discriminate both options, with even a somewhat better correspondence for the turbid medium (Fig. 5, Table 1). Therefore, we will not consider the porous envelop approach in the remaining and we will focus the discussion on the results relative to the turbid medium. Though Bayes factor values are rather similar for “E” and “B” crown types for the sequential optimization (Fig. 3b), BF values show a general decreasing tendency as geometrical complexity in crown representation increases (i.e., from “E” to “M” crown types) for the joint optimization (Fig. 3c). Furthermore, for both optimization strategies, the “M” crown type presents the BF values substantially lower than the other crown types. These observations indicate that the best fits are here obtained for the most complex crown shapes, which contrasts with the optimization using only radiation data for which the best fits occurred for simpler crown forms (i.e., “Ed” and “Bd” crown types) though with rather limited differences in BF values compared with the other crown forms (Fig. 3a). These results demonstrate that accounting for crown asymmetry is much more important to estimate light interception by individual trees than to predict transmitted radiation. Indeed, the radiation reaching a certain point on the ground comes from beams crossing many trees and errors made for some trees are compensated by others. Pleading in favor of the more complex crown shapes, optimizations involving growth data highlight the role of crown plasticity in the aboveground competitive neighborhood interactions and in improving the efficiency of light interception (Longuetaud et al. 2013; Pretzsch 2014; Kru̇ček et al. 2019). It is worth noting that, in the present study, tree growth data were considered as an indicator of the radiation intercepted by tree crowns though, besides radiation, photosynthesis is also influenced by water availability and nutrient status (Walker et al. 2014). Water availability was accounted for in the determination of GPP through its influence on stomatal conductance (see Section 2.1.4), while the effect of nutrient status on carbon assimilation was not considered in the model version used for this study. Nevertheless, comparable foliar nutrient concentrations are observed for the four sites for which growth data were used (data not shown). Therefore, the nutrient status should not have led to contrasted photosynthesis efficiencies among study sites.

The values retrieved for the extinction coefficient k (see Fig. 6) are either lower than (oak), close to (beech) or much higher than (hornbeam) the value of 0.5 assuming spherical distribution of the leaf orientation and widely adopted in 3D crown models (e.g., Brunner 1998; Piboule 2001; Courbaud et al. 2003). Discrepancies between our k estimates and this theoretical value might arise from divergences in real leaf orientation distribution from the theoretical one and/or from clumping of leaves within the crowns, the clumping factor Ω having been set to 1 (see Section 2.3.1). Regarding canopy openness p, our estimates are in the upper range of values found in literature for broadleaved species, ranging from 0.05 to 0.56 (Canham et al. 1994, 1999; Beaudet et al. 2002, 2011; Astrup and Larson 2006; Lefrançois et al. 2008; Da Silva et al. 2012). p is acknowledged to differ according to species, tree size and site characteristics (Canham et al. 1999; Astrup and Larson 2006; Lefrançois et al. 2008) while its estimated values may also depend on the method adopted for their determination, generally either from photographs of isolated crowns (Canham et al. 1999; Lefrançois et al. 2008; Boivin et al. 2011; Da Silva et al. 2012) or from model inversion (Groot 2004). Values retrieved for the NPP to GPP ratio r_NppGpp (Fig. 5) are within the 0.22–0.79 range reported in a recent review on this topic (Collalti and Prentice 2019) and are in agreement with values more specifically mentioned for the species considered in the present study (Kutsch et al. 2005; Nagy et al. 2006; Granier et al. 2008; Wu et al. 2013; Campioli et al. 2015).

5 Conclusions

The calibration exercise performed in this study on a 3D-crown forest radiative transfer model (SAMSARALIGHT) showed contrasted model performances depending on the adopted model configuration options (i.e., crown shape and radiative transfer formulation) and on the optimization strategy (i.e., sequential vs joint optimization) used to combine radiation and growth data information content. A first important outcome is the generally better quality of the fits between observed and predicted values obtained for the turbid medium formulation of the radiative transfer compared with the porous envelop approach, which was especially pronounced for the fit with radiation data. This result arises from the formulation adopted for the porous envelope, which does not account for the beam path within the tree crowns and from the intra-specific variations in leaf area density. For the turbid medium approach, only minor differences in model fit quality appeared among crown shapes when the optimization was carried out on radiation data only, while a general improvement of model performances occurred as crown shape complexity increased when involving the growth data. These findings demonstrate that a few hemispherical photos on the ground are not sufficient to calibrate correctly light interception by trees. They also highlight the role of crown plastivity to optimize light interception and maximize tree growth. Adding growth data is therefore necessary to improve the calibration, especially when the model is then used to predict tree growth in heterogeneous forests.

Appendix

The LAD_model option accounting for the variation of SLA with the vertical position within the canopy when determining LAD is implemented as follows. Two sets of n trees each are selected in the vicinity of the tree for which SLA is to be calculated, the first set gathering the neighboring trees presenting the smallest Hcb while the second set contains the trees with the largest Htot. The number of trees n is defined so that up to N trees per hectare of neighborhood area are selected, it is expressed as a function of the neighborhood radius R (m):

$$ n = \mathrm{round \ up}\left( \frac{N \ {\Pi} \ R^{2}}{10 \ 000}\right) $$

(A.1)

The values assigned to N and R are to be defined by the model user. The average height levels of the bottom of the canopy ($\overline {Hc}_{bottom}$, m) and of the top of the canopy ($\overline {Hc}_{top}$, m) over the tree neighborhood are determined by (see Fig. 8a):

• from the set of trees with the smallest Hcb

  $$ \overline{Hc}_{bottom} = \frac{\sum\limits_{i=1}^{n}(Hcb_{i}+Hlce_{i})/2}{n} $$

  (A.2)

• from the set of trees with the largest Htot

  $$ \overline{Hc}_{top} = \frac{\sum\limits_{i=1}^{n}(Hlce_{i}+Htot_{i})/2}{n} $$

  (A.3)

Fig. 8

Schematized determination of crown SLA for the LAD_model option. (a) $\overline {Hc}_{bottom}$ and $\overline {Hc}_{top}$ are calculated over the tree neighborhood (with n set to 3 here) and (b) crown SLA is then estimated through linear interpolation of SLA_bottom and SLA_top according to crown height level between $\overline {Hc}_{bottom}$ and $\overline {Hc}_{top}$

For crowns without lower and upper parts subdivision (i.e., “Ec” and “Ed” types, and all “P” and “C” crowns), the crown SLA value is obtained by linear interpolation of SLA_bottom and SLA_top according to crown height level using (see Fig. 8b):

$$ SLA = \frac{SLA_{top}-SLA_{bottom}}{\overline{Hc}_{top}-\overline{Hc}_{bottom}} (\overline{Hc}-\overline{Hc}_{bottom}) + SLA_{bottom} $$

(A.4)

where $\overline {Hc}$ (m) is the tree average crown height evaluated as :

$$ \overline{Hc} = (Hcb+Htot)/2 $$

(A.5)

For vertically subdivided crown types (i.e., “Bc,” “Bd” and “M”), crown SLA is determined separately for the lower (SLA_down, m² kg^− 1) and the upper (SLA_up, m² kg^− 1) parts by applying Eq. (A.4) for the average lower crown height level ($\overline {Hc}_{down}$, m) and for the average upper crown height level ($\overline {Hc}_{up}$, m), respectively. $\overline {Hc}_{down}$ and $\overline {Hc}_{up}$ are computed as :

$$ \overline{Hc}_{down} = (Hcb+Hlce)/2 $$

(A.6)

$$ \overline{Hc}_{up} = (Hlce+Htot)/2 $$

(A.7)

These interpolated SLA values are then converted to LAD using estimates of tree leaf biomass (LeafBiomOM, kg of organic matter) provided by an allometric relationship :

$$ LeafBiomOM = \alpha_{leaf} \ C130^{\beta_{leaf}} \ Dd^{\gamma_{leaf}} $$

(A.8)

where Dd (-) is the ratio between the diameter of the crown projection (D, m) determined from the tree mean crown radius and the diameter of the trunk at breast height (DBH, m), and α_leaf, β_leaf and γ_leaf are model parameters.

For “Ec”, “Ed,” “P” and “C” crowns, SLA and LeafBiomOM estimates are combined to evaluate modeled LAD (LAD_model, m²m^− 3) through :

$$ LAD_{model} = \frac{SLA \ LeafBiomOM}{V} $$

(A.9)

where V (m³) is the crown volume determined by applying the ellipsoid, the paraboloid or the cone volume formula depending on the crown form.

For “Bc”, “Bd” and “M” crowns, the distribution of leaf biomass among the different crown parts has to be accounted for. Therefore, for these crown types, LAD_model is evaluated for the lower and the upper crown parts using :

$$ LAD_{model, \ down} = \frac{SLA_{down} \ LeafBiomOM \ (1-UFLB)}{V_{down}} $$

(A.10)

and

$$ LAD_{model, \ up} = \frac{SLA_{up} \ LeafBiomOM \ UFLB}{V_{up}} $$

(A.11)

where V_down (m³) and V_up (m³) are the volumes of the lower and the upper crown parts, respectively. UFLB (-) is a parameter representing the proportion of tree leaf biomass located in the upper crown fraction, it is bounded upwards by the value of the ratio V_up/(V_down + V_up).