Analyzing the Vertical Distribution of Crown Material in Mixed Stand Composed of Two Temperate Tree Species "2279

The material distribution inside tree crowns is difficult to quantify even though it is an important variable in forest management and ecology. The vertical distribution of a relative density index (i.e., vertical profile) of the total, woody, and leafy material at the crown scale were estimated from terrestrial laser scanner (TLS) data on two species, sugar maple (Acer saccharum Marsh.) and balsam fir (Abies Balsamea Mill.). An algorithm based on a geometrical approach readily available in the Computree open source platform was used. Beta distributions were then fitted to the vertical profiles and compared to each other. Total and leafy profiles had similar shapes, while woody profiles were different. Thus, the total vertical distribution could be a good proxy for the leaf distribution in the crown. Sugar maple and balsam fir had top heavy and bottom heavy distributions respectively, which can be explained by their respective architectural development. Moreover, the foliage distribution of sugar maples shifted towards the crown base when it was found in mixed stands, when compared to pure stands. The opposite behavior was observed for balsam firs, but less pronounced. According to the shape of the foliage distribution, sugar maple takes advantages from mixture contrarily to balsam fir. From a methodological point of view, we proposed an original approach to separate wood from leaf returns in TLS data while taking into account occlusion. Wood and leaf separation and occlusion problems are two challenging issues for most TLS-based studies in forest ecology.


Introduction
Quantifying vertical distribution of foliage is of major importance to estimate light interception at the tree crown scale [1,2]. In comparison to the so-called big-leaf approach where the entire tree leaf area operates as a single leaf [3], an explicit vertical distribution of foliage allows a more accurate quantification of the amount of light absorbed and the subsequent tree growth [4][5][6]. The linear relationship between light absorbed and tree growth is now generally accepted [4,5] and is used in several process-based models to understand and predict tree and stand growth. It is also used in forest management for growth and yield models [7][8][9].
These approaches using voxel and allowing to correct occlusions biases have also been used to quantify plant area density (PAD) [24,45,46].
The main objective of our study was to analyze the variability of the material distributions inside the crown between a coniferous (balsam fir: Abies Balsamea Mill.) and a broadleaved (sugar maple: Acer saccharum Marsh.) species in pure and mixed stands at two developmental stages, which are both shade tolerant. We made the hypotheses that (i) the distribution of the balsam fir will be bottom heavy when compared to sugar maple since their intrinsic architectural development (i.e., genetic determinism) differs [47]. Balsam fir, similar to many coniferous species, has a strong apical dominance with horizontal branches compared to sugar maple that loses its apical dominance fairly quickly and produces large forks; (ii) the vertical distribution of the plant material of trees in pure stands will be top heavy when compared to trees in mixed stands due to competition reduction in mixed stands as consequence of crown complementarity [21,22]; and (iii) the distribution of mature trees will be top-heavy for sugar maples and bottom heavy for balsam firs when compared to younger trees, due to their differences in crown architecture, i.e., higher insertion angles of the branches for sugar maple, when compared to the nearly horizontal angle of balsam fir branches. We tested the hypothesis that (iv) simple crown attributes explain the differences in vertical distributions between species, stand type and developmental stage. In order to answer this main objective, an approach to estimate the vertical distribution of a tree's material using TLS data was added as a methodological objective. More specifically, we propose to (1) apply a geometric approach to separate woody and leafy returns in the TLS point cloud and (2) apply a method to deal with the limitations of signal occlusion when using TLS data in forest.

Sites and Tree Sampling
A total of 36 balsam fir and 36 sugar maple target trees were sampled among six sites located in eastern Quebec, Canada, as described in Martin-Ducup et al. [48]. Provincial eco-forest maps [49] were used to select the sites, based on homogenous edaphic and slope characteristics. Each site was composed of three stand types: a pure balsam fir stand, a pure sugar maple stand, and a mixed stand where both species were found. In each site, three codominant balsam fir trees were randomly selected in the pure balsam fir stand and three in the mixed stand (i.e., 3 balsam fir trees in pure stand and 3 balsam fir trees in mixed stand per site X 6 sites = 36 balsam fir trees). In the same way, three codominant sugar maple trees were randomly selected in the pure sugar maple stand and three in the mixed stand (i.e., 3 sugar maple trees in pure stand and 3 sugar maple trees in mixed stand per site X 6 sites = 36 sugar maple trees). Sites differed from each other by the developmental stage of the co-dominant trees and can be divided into two groups: intermediate and mature stage sites. Tree characteristics by site and by species are given in Table 1. Basal areas were estimated using four wedge prism plots per stand type. Mean basal area by species for each condition and its developmental stage are given in Table 2.  Target trees (TT) were scanned with a TLS, the FARO ® Focus 3D (FARO Technologies Inc., Lake Mary, FL, USA), during the summer of 2013 at four opposite viewpoints after the neighbor trees were felled. The TTs were scanned with a TLS with a resolution angle of 0.036 • in both the horizontal and vertical directions. The scanner rotation angles ranged from 0 • to 360 • in the horizontal direction and from 60 • to 90 • in the vertical direction. The FARO Focus 3D TM is a single return instrument. Scan alignment was required to create one global point cloud from the four scans. Alignment was facilitated by using six spherical targets placed around the TT, with a minimum of three targets visible in each scan. No filters were activated on the FARO instrument during the scans, but two filters were applied during data preprocessing using the software FARO Scene 5.0 (FARO Technologies Inc., Lake Mary, FL, USA): (i) the dark scan points filter to remove points with a very low reflectance (less than 300) which correspond to false returns or sky points and (ii) the stray filter, which validates that the return is due to an object by comparing the distance between the object and the scanner with the distance of the surrounding points (grid size of 3 × 3 × 3 points). The return is retained as a scan point if 50% of the differences in distance with the surrounding points is less than 0.02 m.
Files were prepared for wood/leaf separation and the correction of signal occlusion by extracting several files format. An ASCII file containing the four aligned and filtered point clouds was exported from FARO ® Scene to be used by the wood/leaf discrimination algorithm. The TT of each site was previously manually isolated using the Computree software (http://rdinnovation.onf.fr/, consulted on 14 March 2017). Tree isolation was easy since the TT neighbors had been felled before the scans were taken and the TTs were thus left clearly isolated in the middle of a gap (Figure 1a). Filtered point clouds were necessary for this step since raw point cloud contained too much noise points which could bias the geometrical approach. A raw point cloud file ".xyb" without filters was extracted for each scanning position to be used by the occlusion correction algorithm (see section Voxel representation of material distribution). The ".xyb" file was chosen as it is faster to extract from FARO ® Scene and faster to process with the software used, Computree (http://rdinnovation.onf.fr/, consulted on 14 March 2017). Target trees (TT) were scanned with a TLS, the FARO ® Focus 3D (FARO Technologies Inc., Lake Mary, FL, USA), during the summer of 2013 at four opposite viewpoints after the neighbor trees were felled. The TTs were scanned with a TLS with a resolution angle of 0.036° in both the horizontal and vertical directions. The scanner rotation angles ranged from 0° to 360° in the horizontal direction and from 60° to 90° in the vertical direction. The FARO Focus 3D TM is a single return instrument. Scan alignment was required to create one global point cloud from the four scans. Alignment was facilitated by using six spherical targets placed around the TT, with a minimum of three targets visible in each scan. No filters were activated on the FARO instrument during the scans, but two filters were applied during data preprocessing using the software FARO Scene 5.0 (FARO Technologies Inc., Lake Mary, FL, USA): (i) the dark scan points filter to remove points with a very low reflectance (less than 300) which correspond to false returns or sky points and (ii) the stray filter, which validates that the return is due to an object by comparing the distance between the object and the scanner with the distance of the surrounding points (grid size of 3 × 3 × 3 points). The return is retained as a scan point if 50% of the differences in distance with the surrounding points is less than 0.02 m.
Files were prepared for wood/leaf separation and the correction of signal occlusion by extracting several files format. An ASCII file containing the four aligned and filtered point clouds was exported from FARO ® Scene to be used by the wood/leaf discrimination algorithm. The TT of each site was previously manually isolated using the Computree software (http://rdinnovation.onf.fr/, consulted on 14 March 2017). Tree isolation was easy since the TT neighbors had been felled before the scans were taken and the TTs were thus left clearly isolated in the middle of a gap ( Figure 1a). Filtered point clouds were necessary for this step since raw point cloud contained too much noise points which could bias the geometrical approach. A raw point cloud file ".xyb" without filters was extracted for each scanning position to be used by the occlusion correction algorithm (see section Voxel representation of material distribution). The ".xyb" file was chosen as it is faster to extract from FARO ® Scene and faster to process with the software used, Computree (http://rdinnovation.onf.fr/, consulted on 14 March 2017).

Wood and Leaf Separation
A cylinder fitting algorithm was used to determine which returns were tied to a wood structure and which ones were from a leaf. This algorithm is available in the Computree software and was initially designed to detect stems. In the present study we kept the default parameters proposed in the software. The approach also needs some treatments in R (cylinder filtering in step 4 step 5) [49]. The algorithm processes each point cloud to identify the linear and cylinder-like features by applying the following five steps: (1) Clustering of the point cloud into 2 cm high horizontal slices in order to identify small trunk arcs and branches. Two points are considered part of the same cluster if they are separated by no more than 3 cm.
(2) Vertically merging point clusters to form linear elements following these two conditions: (i) the vertical distance between clusters has to be less than 10 cm and (ii) the bounding box of two clusters must intersect each other in the XY plan. This step was originally designed to detect tree stems, but a continuous branch which is not vertical may also meet the two conditions.
(3) Merging parallel linear elements to form the trunk and branches. The XY distance threshold to merge linear elements must be large enough to merge non-vertical branches. It was set to 50 cm after preliminary visual tests on the datasets.
(4) Adjusting cylinders on linear elements by using variable angles (from vertical to horizontal). This adapts the cylinder fitting procedure on tree stem segments and branches that are not vertical. Several cylinders can be adjusted for each linear element. The points associated to a cylinder are projected on a plan, perpendicular to the axis of the cylinder, and a circle is fitted on the 2D point cloud. An absolute and relative error of the fitted circle is associated to each cylinder. Moreover, a linearity error is also associated to each cylinder. This last parameter is the mean square error of the 3D linear regression on the barycenter of three consecutive cylinders (i.e., belonging to the same linear element). Cylinders considered as a wood linear element in the point cloud has a maximum absolute and relative circle error of 4 cm and 30%, respectively, and a linearity error of 0.001 m.
More details on each step are given in the Computree software documentation (http:// rdinnovation.onf.fr/, consulted on 14 March 2017).
(5) Exporting the point cloud with an associated cylinder ID. The points associated with a cylinder (points of step 4) were considered as returns from wood (R W ) and the other points were considered to be returns from leaves (R L ). Finally, the point cloud was divided into 10 cm edged voxels to calculate the proportion of R W and R L in each voxel ( Figure 1).

Voxel Representation of Material Distribution
The number of returns in the point cloud of a forest scene is not an exact representation of the 3D distribution of the canopy components because TLS data is biased due to signal occlusion and variable sampling rates. An algorithm called L-vox, available in Computree, allows matching more than 95% of the distributed material between simulated scans on a simulated tree and the actual distribution of the simulated tree's material [42]. The general idea of the L-vox algorithm is, firstly, to divide the point cloud into a matrix composed of voxels (cubes with edges of equal size). Secondly, to calculate a total Relative Density Index (RDIT) for each voxel using the number of laser beams that crosses the voxel and the number of returns from within the voxel. Total Relative Density Index appellation was used because no differentiation between woody and leafy material was done at this stage. The algorithm assumes that the TLS is a single return instrument. The theoretical beams are calculated using a ray-tracing model based on the scan parameters (angular resolution, scanner position, and vertical and horizontal openness angle). The RDIT is thus calculated as follows: where n r is the number of observed returns in the voxel, n t the potential total number of beams passing through the voxel (if no objects existed between the voxel and the scanner) obtained from the ray-tracing algorithm, and n b the number of intercepted beams between the scanner and the voxel estimated using the TLS data and the ray tracing algorithm (see [50] for further explanations and a schematic representation). The algorithm was applied to each of the four scans, but the four voxel-grids of each target tree (TT) were combined for analysis. This implies that each voxel has four RDIT values, each corresponding to a scan. The retained RDIT value is the one corresponding to the scan that has the least occlusion (i.e., the highest n t − n b value). The voxel edge size was 10 cm, as recommended in Béland et al. [24] and Fournier et al. [42], corresponding to a voxel volume (VoxV) of 0.001 m 3 .
In order to extract all the voxels belonging to the TT, the isolated TT point cloud previously described was used to select the voxels from the L-vox grid. More precisely, the crown boundaries of the isolated TT point cloud were identified using 2D convex hull for every 10 cm high slice with the "Geometry" R package [51]. The resulting polygons were used to identify all the voxels of the L-vox grid inside those polygons. These voxels were considered to belong to the TT.
At this stage, two TT voxel grids are available: (1) the R W /R L voxel grid where each voxel has the proportion values of R W (and thus of R L ), and (2) the L-vox grid, where each voxel has its RDIT value. The RDIT of the voxel in the L-vox grid was then divided into wood RDI (RDIW) and leaf RDI (RDIL). This step was done by crossing the wood/leaf proportion voxel grid with the RDIT voxel grid using the following equations ( Figure 2).
The resulting three voxel grids (i.e., RDIT, RDIW, and RDIL) were used to draw vertical profiles. Since the aim of the study is to analyze crown material distribution, only the voxels of the crown (the tree stem below the crown base was excluded) were used for the three profiles. Cumulative relative RDI at each vertical layer (∆z = 10 cm) were used to make the profiles (RDI in Equations (4) and (5) can be replaced by RDIT, RDIW, or RDIL): where RDI z is the RDI value in the vertical layer z, nz is the number of voxel in the layer of height Z and nt the number of voxel of the entire tree crown. The shape of the profile was modelled using a two parameters (r, s) beta distribution. This distribution was chosen because it is highly flexible and defined within an interval of fixed endpoints [10,12]. The beta distribution (Equation (3)) was fitted to each tree and density distribution (RDIT, RDIW, RDIL) individually using the nls function in the R software [52]: The RDIT voxel grid was also used to quantify the following six crown metrics ( Figure 3): • Crown length (CrL): The highest Z coordinate was used to determine the height of the tree. The Z position of the crown base was recorded manually (in FARO ® Scene) by using the aligned point clouds for each tree. • Projected area of the crown (PACr): Area of the polygon calculated on the XY plane using the 2D convex hull algorithm of the "Geometry" package [51]. • Crown volume (CrV): Volume calculated from the hull generated by the 3D alpha shape with the R package "Alphashape3d" [53].

•
Crown density (CrD): Proportion of crown volume occupied by material: • Crown width to length ratio (CrW:CrL): The crown width is estimated from the diameter of a circle whose area is equal to the projected crown area.

•
Crown proportion (Cr%): Crown length to tree height ratio. The RDIT voxel grid was also used to quantify the following six crown metrics ( Figure 3): • Crown length (CrL): The highest Z coordinate was used to determine the height of the tree. The Z position of the crown base was recorded manually (in FARO ® Scene) by using the aligned point clouds for each tree. • Projected area of the crown (PACr): Area of the polygon calculated on the XY plane using the 2D convex hull algorithm of the "Geometry" package [51]. • Crown volume (CrV): Volume calculated from the hull generated by the 3D alpha shape with the R package "Alphashape3d" [53]. • Crown density (CrD): Proportion of crown volume occupied by material: • Crown width to length ratio (CrW:CrL): The crown width is estimated from the diameter of a circle whose area is equal to the projected crown area. • Crown proportion (Cr%): Crown length to tree height ratio.

Data Analysis
Data analysis was divided into two parts. First, we analyzed the effect of species, stand type, and developmental stage on the r and s parameters of the beta distribution using mixed effect models. Second, we modelled the linear relationships between each parameter of the beta distribution and the crown metrics, as well as between the parameters and the stand characteristics. These two models were compared.
The effects of stand type, developmental stage and species on the r and s parameters were analyzed using the following model for the three types of profile (i.e., RDIT, RDIW, and RDIL): where is the r or s beta distribution parameter (of RDIT, RDIW or RDIL) estimated for a species (k) of stand type ( ) in the site ( ), and the overall mean. The variable is a binary variable indicating stand type ( = 0 for mixed stand, = 1 for pure stand), is a binary variable that indicates the developmental stage ( = 0 for intermediate stage and = 1 for mature stage), is an indicator of the species ( = 0 for sugar maple and = 1 for balsam fir), is the site random effects ( ~ (0, )) and the residual error ( ~ (0, )).
All possible combinations of presence/absence of the crown metrics and stand characteristics (i.e., SC that could be SP, ST, or SP) were then calibrated (Equation (6)), and a total of 479 models

Data Analysis
Data analysis was divided into two parts. First, we analyzed the effect of species, stand type, and developmental stage on the r and s parameters of the beta distribution using mixed effect models. Second, we modelled the linear relationships between each parameter of the beta distribution and the crown metrics, as well as between the parameters and the stand characteristics. These two models were compared.
The effects of stand type, developmental stage and species on the r and s parameters were analyzed using the following model for the three types of profile (i.e., RDIT, RDIW, and RDIL): where Y ijk is the r or s beta distribution parameter (of RDIT, RDIW or RDIL) estimated for a species (k) of stand type (j) in the site (i), and µ the overall mean. The variable ST ij is a binary variable indicating stand type (ST ij = 0 for mixed stand, ST ij = 1 for pure stand), DS ij is a binary variable that indicates the developmental stage (DS i = 0 for intermediate stage and DS i = 1 for mature stage), SP ijk is an indicator of the species (SP ijk = 0 for sugar maple and SP ijk = 1 for balsam fir), u i is the site random effects ( u i ∼ N 0, σ 2 i ) and ε the residual error ( ε ijk ∼ N 0, σ 2 i ).
All possible combinations of presence/absence of the crown metrics and stand characteristics (i.e., SC that could be SP, ST, or SP) were then calibrated (Equation (6)), and a total of 479 models compared. Among the six metrics that were initially calculated, crown projected area (PACr) was dropped because it was found to be collinear with other variables (especially with CrV), as indicated by a high variance inflation factor (VIF > 5) [54]. Once the PACr dropped, the VIF values of the remaining variables were below five. Interaction effects between the stand characteristics were also considered when two of the three stand characteristics (SC) were in the same model: where a x are the fixed effect parameters for the crown metric x (CM x ); SC x is a stand characteristic variable x among the three characteristics tested in the first step (ST, DS or SP), u i is the site random effects parameters (u i ∼ N 0, σ 2 i ).

Species, Type of Stand and Developmental Stages Profile Variations
The total material vertical distribution for balsam fir was bottom heavy, whereas the opposite was observed for sugar maple. This was translated by a mode around 75% of the height within the crown for sugar maple and around 30% for balsam fir (Table 3, Figures 4 and 5). 3D representation and associated RDIT profiles of four mature example trees are illustrated in Figure 6. compared. Among the six metrics that were initially calculated, crown projected area (PACr) was dropped because it was found to be collinear with other variables (especially with CrV), as indicated by a high variance inflation factor (VIF > 5) [54]. Once the PACr dropped, the VIF values of the remaining variables were below five. Interaction effects between the stand characteristics were also considered when two of the three stand characteristics (SC) were in the same model: where are the fixed effect parameters for the crown metric ( ); is a stand characteristic variable x among the three characteristics tested in the first step (ST, DS or SP), is the site random effects parameters ( ~ (0, )).

Species, Type of Stand and Developmental Stages Profile Variations
The total material vertical distribution for balsam fir was bottom heavy, whereas the opposite was observed for sugar maple. This was translated by a mode around 75% of the height within the crown for sugar maple and around 30% for balsam fir (Table 3, Figures 4 and 5). 3D representation and associated RDIT profiles of four mature example trees are illustrated in Figure 6.    On the upper part, the voxel of the crown is represented with a color gradient from yellow to red for RDIT values (similar to Figure 2). The grey parts of the trees are the trunk and the dead branches (not considered for the profile). On the lower part the RDIT profile of each tree is represented.   On the upper part, the voxel of the crown is represented with a color gradient from yellow to red for RDIT values (similar to Figure 2). The grey parts of the trees are the trunk and the dead branches (not considered for the profile). On the lower part the RDIT profile of each tree is represented. On the upper part, the voxel of the crown is represented with a color gradient from yellow to red for RDIT values (similar to Figure 2). The grey parts of the trees are the trunk and the dead branches (not considered for the profile). On the lower part the RDIT profile of each tree is represented. Table 3. Analysis of variance for both beta distribution parameters (r and s) and the three profile types (Equation (5)). Fixed effect of stand type (ST), developmental stage (DS) species (SP), the interaction between stand type and developmental stage (ST:DS), between stand type and species (ST:SP) and between developmental stage and species (DS:SP) are shown. In mixed stands, the vertical distribution for balsam fir was top heavy relative to pure stand, while the opposite was observed for sugar maple. The effect of stand type was significant on s parameter with a strong interaction between species and stand type on both parameters (Table 3). For mature sugar maple, in mixed stands the mode was located closer to the base of the crown and the variance of the distribution was higher relative to pure stands (Figure 4). For balsam fir in mixed stands, the mode was located higher in the crown relative to pure stands for both the mature and intermediate developmental stages (Figure 5).

RDI
Mature balsam fir trees had a distribution of material lower in the crown than younger trees contrarily to sugar maple (Figures 4 and 5). The effect of developmental stage was significant for the s parameter as well as the interaction between species and developmental stage on the r parameter for RDIT and r and s parameters for RDIL (Table 3).

Crown Metrics Models
Model comparisons generally showed that crown metrics were not sufficient to explain the variability of the distribution profiles. Nevertheless, two crown metrics (i.e., crown density and crown proportion) were retained in all the best models ( Table 4). The best models for both r and s parameters were the same for every distribution considered (i.e., RDIT, RDIW or RDIL). Crown density (CrD) was only significant for the s parameters of RDIT and RDIL distribution (p < 0.0001) although crown proportion (Cr%) was significant for all parameters of all distributions (p ranged between <0.0001 and <0.026) except for s parameters of the RDIW distribution (p = 0.46). All stand characteristics (species, stand type and developmental stage) as well as their interactions were retained for the r parameters. All interaction terms were found to be statistically significant for RDIT and RDIL (p ranged between <0.0001 and <0.0059). For RDIW, only the interactions between stand type and species (p = 0.03) and stand type and developmental stage (p = 0.019) were statistically significant. For s, the only stand characteristic that remained in the models was the species, which was significant in all the models (p ranged between <0.0001 and <0.0012). Table 4. Estimates and standard error for the best models to predict r and s beta distributions parameters (Equation (6)). AIC criteria were used to select each best model among 479.

Type of Material Comparison
Total material profile and leaf profile were similar whereas wood profile was lower for both species. For the sugar maple, the RDIW distribution had a higher variance and the mode was located lower in the crown than for RDIT or RDIL. For balsam fir, the RDIW mode was higher in the crown with a higher variance than for the RDIT or RDIL distributions. The only significant effect on RDIW parameters, other than differences between the species, was the interaction between stand type and species for the parameter r ( Table 3). The effect of stand type is stronger for sugar maple than for balsam fir, with a lower mode in mixed stands than in pure stands (Figures 4 and 5). This effect was more pronounced for mature trees although there was no significant effect of the interaction between the developmental stage and stand type (Table 3).

Differences between Species
Differences between the two species under study were found, and these can be explained by their crown architecture. Sugar maple trees had top-heavy distributions and balsam fir trees bottom heavy ones, for every kind of return considered (total, wood, or leaves) (Figures 4 and 5). Balsam fir conforms to the Massart architectural model, where a strong apical dominance with a well-defined unique trunk and horizontal branches are observed [47]. This architectural development is constant across the ontogeny and leads to a conical, bottom heavy crown (Figures 6 and A1 in Appendix A). Sugar maple, however, generally has multiple forks due to the early loss of apical dominance in its ontogeny. Co-dominant trees generally present an inverted cone shape with a high proportion of material towards the top of the crown [55] (Figures 6 and A2). Moreover, the difference in the architectural development of these two species also explains the observed opposite response of the developmental stage. The high insertion angle on the trunk (i.e., horizontal) of the balsam fir's first order branches leads to an increase in material low in the crown. On the other hand, the low insertion angle (i.e., vertical) of the sugar maple's first order branches leads to an accumulation of material higher in the crown (Figures A1 and A2). This became more pronounced in older sugar maples.
Coniferous species usually concentrate their foliage in the lower part of the crown, although some inter-species variations has been explained by species' shade tolerance. Our results on balsam fir agreed with previous studies, where a strong concentration of foliage around 40% of the crown height was observed [11,13] (Figure 5). On the other hand, studies on light demanding pine species highlighted that jack pine (Pinus banksiana Lamb.) has a concentration of foliage in the upper part of the crown [56] and scots pine (Pinus sylvestris L.) close to the midpoint [57].
Few studies were conducted on the vertical distribution of foliage in broadleaved species [58], and even fewer compare broadleaved and coniferous. As with conifers, the foliage distribution of broadleaved species is strongly correlated to shade tolerance. The foliage distribution of shade tolerant species is bottom heavy compared to shade intolerant [15,16,59]. Considering that the sugar maple is very tolerant to shade, it was surprising to observe such a top-heavy material distribution (Figure 4). An explanation could be the fact that in this study, we defined the crown base as the insertion height of the first fork or main branch that has leaves (Figures 6 and A2). Therefore, the bottom part of the crown is generally composed by two main branches with a few twigs and leaves. Higher in the crown, the branching topology is more complex, with more twigs and leaves and thus, more material sampled by the TLS. These results were confirmed by the differences in the shapes of wood distribution when compared to total (or leaf) distribution ( Figure 4). The distribution peak for wood returns was located lower in the crown than for the total (or leaf) returns. This is consistent with the hypothesis that big branches with very few twigs and leaves mostly occupy the bottom of the crown. Finally, these differences between wood and total (or leaf) distributions were more important for mature trees than for intermediate ones. This is in accordance with the fact that sugar maple trees have generally smaller branches low in the crown at young developmental stages [47].

Differences between Mixed and Pure Stands
Sugar maple had a more homogeneous distribution lower in the crown in mixed stands when compared with sugar maple trees in pure stands (Figure 4). A downward shift of the material is often associated to trees that are confronted to less competition. It is advantageous in terms of light interception because this minimizes self-shading, thus increasing light interception [10,11,14]. This is in accordance with the results of Martin-Ducup et al. [21], which were obtained from the same sites, where they showed that the competitive pressure of sugar maples is lower and space occupancy is higher in mixed stands than in pure ones. Metz et al. [60] also showed that the competition is lower in mixed stands, which in turn explains the higher growth of European beech.
It can be difficult to separate the effect of diversity from the effect of competition in natural stands. However, it was interesting to observe that in the case of sugar maple, stand density (i.e., competition) is higher in mixed stands than in pure ones ( Table 2). The advantage of mixture seems to overcome the potential negative effect of density (i.e., competition). Garber and Maguire [10] observed that the positive effect of diversity on vertical distribution was more important in high density stands. As stand density increases, the interaction between the different individuals is stronger. This increased interaction thus underlines the importance of trait complementarity in mixed stands [61][62][63]. The foliage of balsam firs was concentrated higher in the crown in mixed stands than in pure ones ( Figure 5). This suggests that in mixed stands, they experience more competition and have less efficient crowns to intercept light. These results were true for both mature and intermediate developmental stages. Here again these differences could not be attributed to stand density because the basal area of the balsam fir was higher in pure stand than in mixed stand, contrarily to the sugar maple (Table 2). By considering the stem density only, one would have expected an inverse response. But here, the stand composition seems to be in cause.
When balsam fir and sugar maple were found in the same stand, it appears that one species takes advantage of the mixture to the detriment of the other. Both species had opposite reactions to mixture, with balsam firs shifting their foliage up towards the middle of the crown and sugar maples shifting it down towards the bottom of the crown. Mixture seems to be advantageous for sugar maples but not for balsam firs. However, the distribution was generally more homogenous along the crown height (i.e., higher variance) in mixed stands for both species and could result in a better total space occupation with an optimal sharing of the available space in mixed stands. Visual observation on the graphical representations of balsam fir trees ( Figure A1) suggest that crowns occupy more space in mixed stand (at least for mature trees). An analyzes of the global space occupation in response to mixture, such as in Martin-Ducup et al. [21] on sugar maple, can provide some further support for these observations. Finally, our stands are composed of many other species. An analysis of the response of all species is needed to clarify the effect of mixture on space occupation at the stand scale.

The Effect of Developmental Stage
The effect of stand composition observed for sugar maple was only significant for mature trees, as presented in both Martin-Ducup et al. [21] and the present study ( Table 3). The potential benefit of increased diversity varies with the developmental stage due to temporal modifications in resource capture and utilization and this could partly explain our results [15,64]. Moreover, one could argue that the stronger effect of mixture on mature trees was confounded with the effect of stand density. However, the density of intermediate and mature trees in mixed stands were very close ( Table 2).

Do Crown Metrics Explain Material Distribution?
The overall shape and space occupation of the crown can be good proxies of the leaf area and leaf distribution in the crown [14,55]. We quantified the effect of several simple and complex crown metrics on the profiles shape. CrD and Cr% were the crown variables that best explained the variability of the beta distribution parameters. The Cr% was present and significant in almost all models (Table 4). This suggests that it might not be necessary to quantify complex metrics from TLS data such as PACr, CrV or CrD to obtain the vertical distribution of material. Although CrD was present in all models, it was often non-significant, or had an F-value that was lower than the one observed for the Cr%. However, contrary to our expectations, the crown metrics were not sufficient to explain the variability of the distribution parameters (Table 4). In other words, the differences in vertical distribution between stand type, species and developmental stage cannot be solely explained by crown metrics such as crown shape, dimension, or density. The local competition for light prior to release might better explain the variability of material distribution, as it was the case for space occupation crown metrics for sugar maples in Martin-Ducup et al. [21].

Total and Leaves Distribution
Total material distribution is a good proxy of leafy material distribution. Indeed, the distribution of leaves and total material were found to be very similar, regardless of species. This result could have important implication considering the uncertainties, the algorithm complexity, the time, and the cost associated with wood and leaf separation. The gain obtained by separating wood and leaves to assess material distribution does not appear important, considering how the total and leaf curves are similar, at least for the stands studied in the present paper (Figures 4 and 5). Moreover, with ALS data, occlusion from above and the lower return density do not allow accurate branch representations [43]. Therefore, it seems that the present geometrical method could be biased for large areas using ALS data. However, if our results are confirmed to be independent on potential misidentification of leaves by other studies for different stands and species, characterization of the vertical material profiles using a normalized distribution of all the ALS returns (e.g., using RDIT values from L-vox) would be a good proxy to estimate the leaf profile on large area.

The Advantages and Limitations of the Method
Wood/leaf separation using geometrical shape (i.e., cylinder fitting) presents an advantage over other methods as it can be applied to a point cloud with xyz coordinates information without requiring information on return intensity or color. Geometrical-based methods are more generalizable (i.e., they can be applied to any sensor) and have already proven to be efficient on adult trees [33,38]. Moreover, in the present study, the cylinder fitting algorithms used are already available and operational to detect stems in Computree, an open software (http://rdinnovation.onf.fr/, consulted on 14 March 2017).
The wood/leaf separation approach is a good starting point using a user-friendly geometrical approach on TLS data. It can however be improved and has to be validated. We did not validate the approach since no leaf data was available for the experimental trees. In a further step, optimization of certain parameters of the algorithm steps for each species should be performed. For example, species-specific thresholds may need to be established for the cylinder fitting steps (e.g., cylinder angle; size; error, etc.). Moreover, the cylinder-fitting algorithm has been developed to detect the trunk and is optimized for more or less vertical structures. It can detect horizontal branches, but less efficiently, and thus could induce a large bias for species such as balsam fir. Adding another procedure to the first step of the algorithm to consider the angular direction of the linear clusters could improve the efficiency of horizontal branch detection.
Signal occlusion and variable sampling density are known biases for TLS data, and lead to two major limitations. The first is the underestimation of the material density. Using the L-vox algorithm, we were able to limit this problem to the point where the RDI value of the voxel was assumed to approach an accurate 3D distribution of material [42]. The second limitation is the detection of cylinders in zones with strong occlusion. Although the L-vox algorithm is able to re-estimate the density in a voxel by correcting the occlusion, the spatial arrangement of points is unpredictable. We supposed that this was not a major problem for sugar maples, but it could induce important biases for balsam firs. Indeed, balsam firs have denser crowns and the branch occlusion could be greater than for the sugar maple. For both species, the most occluded part is on small branches that could be assigned as leaves. However, the impact of these small branches on the vertical profiles should be minimal. In order to completely validate the approach used here, it would be interesting to demonstrate the reliability of the approach by running the algorithms on simulated tree [27,65].
Finally, we used TLS data to quantify the vertical distribution of the relative density of the material but we did not transform the density to an area unit (e.g., leaf area). If tree growth and photosynthesis have to be simulated using leaf surface area, a thorough validation of the algorithm is needed to establish if any bias remains. A transformation from the RDIT or RDIL to the leaf area is possible at the voxel scale using methods such as those developed by Hosoi and Osama [26], Van der Zande et al. [66], or Béland et al. [28]. It is also possible to use the total leaf area index at the tree or stand scale obtained by other methods (e.g., field data, allometric relationship), and redistribute it proportionally to RDIL values in order to obtain leaf area at the voxel scale.

Conclusions
We compared the vertical distribution of material inside the crown between a coniferous and a broadleaved species was carried out using data from two types of stands and two developmental stages. The vertical distribution of sugar maple was substantially different than balsam fir. These differences were in accordance with their architectural development along their ontogeny, and more generally, with the differences between conifers and broadleaved species architecture. Their responses to stand composition were opposite, the material distribution of sugar maple was lower in the crown in mixed than in pure stands, contrary to balsam fir. This suggests that sugar maples are facilitated in mixed stand although balsam fir experience more competition than in pure stand. From a methodological point of view, we proposed a geometrical approach for wood/leaf segregation on TLS data based on cylinder fitting while taking into account the occlusion. Wood leaf separation and occlusion are two important challenges in using TLS data in forest ecology. This approach is thus a starting point for a generalizable method to separate wood and leaves using LiDAR technologies and already operational tools (Computree). For one species, total and leaf material vertical distributions were similar, suggesting that the leaf vertical distribution can be inferred from the total distribution of the canopy material that is easier to quantify at large scale using airborne LiDAR data. Figure A2. Three dimensional representation of all sugar maple trees (72 trees). The voxels of the crowns are represented with a color gradient from yellow to red for RDIT values (similar to Figure 2). Figure A2. Three dimensional representation of all sugar maple trees (72 trees). The voxels of the crowns are represented with a color gradient from yellow to red for RDIT values (similar to Figure 2).