Inﬂuence of 3D Spruce Tree Representation on Accuracy of Airborne and Satellite Forest Reﬂectance Simulated in DART

: Advances in high-performance computer resources and exploitation of high-density terrestrial laser scanning (TLS) data allow for reconstruction of close-to-reality 3D forest scenes for use in canopy radiative transfer models. Consequently, our main objectives were (i) to reconstruct 3D representation of Norway spruce ( Picea abies ) trees by deriving distribution of woody and foliage elements from TLS and ﬁeld structure data and (ii) to use the reconstructed 3D spruce representations for evaluation of the effects of canopy structure on forest reﬂectance simulated in the Discrete Anisotropic Radiative Transfer (DART) model. Data for this study were combined from two spruce research sites located in the mountainous areas of the Czech Republic. The canopy structure effects on simulated top-of-canopy reﬂectance were evaluated for four scenarios (10 × 10 m scenes with 10 trees), ranging from geometrically simple to highly detailed architectures. First scenario A used predeﬁned simple tree crown shapes ﬁlled with a turbid medium with simpliﬁed trunks and branches. Other three scenarios used the reconstructed 3D spruce representations with B detailed needle shoots transformed into turbid medium, C with simpliﬁed shoots retained as facets, and D with detailed needle shoots retained as facets D . For the ﬁrst time, we demonstrated the capability of the DART model to simulate reﬂectance of complex coniferous forest scenes up to the level of a single needle (scenario D ). Simulated bidirectional reﬂectance factors extracted for each scenario were compared with actual airborne hyperspectral and space-borne Sentinel-2 MSI reﬂectance data. Scenario A yielded the largest differences from the remote sensing observations, mainly in the visible and NIR regions, whereas scenarios B , C , and D produced similar results revealing a good agreement with the remote sensing data. When judging the computational requirements for reﬂectance simulations in DART, scenario B can be considered as most operational spruce forest representation, because the transformation of 3D shoots in turbid medium reduces considerably the simulation time and hardware requirements. administration, R.J. and L.H.; Resources, J.H.; Software, R.J., N.L. and J.-P.G.-E.; Supervision, L.H., Z.M. and J.-P.G.-E.; Validation, R.J.; Visualization, R.J.; Writing—original draft, R.J., L.H. and Z.M.; Writing—review & editing, L.H., Z.M., J.H., N.L. and J.-P.G.-E.


Introduction
Forests play a key role in the global environment as hotspots of terrestrial biodiversity [1,2], sinks of atmospheric carbon [3], and providers of multiple ecosystem services for human societies [4,5]. Green foliage is the main interface between a tree and its above-ground environment. It is involved in several key ecosystem processes, such as photosynthesis triggered by light interactions (i.e., absorption and scattering), as well as gas and water exchanges [6]. To understand the process-based functional and structural forest changes is a long-term scientific objective [7], and it can be supported by a multitude of optical remote sensing (RS) data and methods e.g., [8][9][10][11][12][13][14].
Optical RS methods require good understanding of the interactions between incoming solar radiation and the vegetation media that are dominated by photosynthetically active foliage. Canopy structure has a particular influence on canopy reflectance characteristics, and especially in the near and short-wave infrared and microwave optical domains [15][16][17]. Although light-canopy interactions can be studied effectively using physically based leaf and canopy radiative transfer models (RTMs), all RTMs rely on a certain level of simplification of vegetation canopy representations. Models simulating vegetation as a formation of horizontally homogeneous layers, such as the family of SAIL models by Verhoef [18] and Verhoef and Bach [19], the semi-discrete model of Gobron et al. [20], and the flux decomposition model of Kallel [21], are better suited to such structurally homogeneous canopies as those of crops, grasslands, or even tree crowns under the assumption of very high canopy closure and foliage density. Better suited to spatially heterogeneous and structurally complex forest canopies, however, are 3D RTMs working with tree architectures created from geometric shapes, such as the Flight model of North [22], FRT model of Kuusk and Nilson [23], or FLiES model of Kobayashi and Iwabuchi [24]. Steadily increasing computational capabilities allow for radiative transfer simulations of highly complex 3D canopies wherein the level of detail reaches to the individual leaf or needle. Such examples include the Raytran model of Govaerts and Verstraete [25], radiosity-graphics model of Qin and Gerstl [26], scalable plant RT model of Bailey et al. [27] and Discrete Anisotropic Radiative Transfer (DART) model of Gastellu-Etechgorry et al. [28,29].
In this study, we used the DART model, which simulates in-situ and RS acquisitions (LiDAR data and images in the optical and thermal spectral domains) and also the radiative budget of any landscape configuration, including terrain topography and atmosphere above the simulated scene [29]. Currently, DART enables constructing representations of trees by (i) creating simple predefined geometric shapes internally, (ii) importing external user-defined 3D objects, or (iii) transforming the imported 3D objects into simplified geometrical shapes [30]. The first of these utilities DART's intrinsic functionalities defining a tree trunk and branches as simple 3D representations. Tree crown is represented by geometric shapes (e.g., an ellipsoid or a truncated cone) that are filled with triangles representing leaf facets or are created by cells known as voxels (i.e., 3D pixels) filled with a turbid medium representing foliage (Figure 1a). The turbid medium represents an infinite number of infinitely small leaves. The triangular facets and the turbid medium are characterized by predefined leaf optical and structural properties. Structural properties are primarily controlled by leaf area index (LAI) and leaf angle distribution (LAD), both specified either per tree or per canopy of the entire simulated scene. Other controlling parameters are distribution of empty voxels (i.e., air gaps) and leaf volume density inside a voxel, which can be specified either per crown or according to crown vertical level. Such a tree representation is computationally very effective and well suited for forest canopies with relatively simple tree architecture, such as those of broadleaf trees. This parametrization may be less suitable for conifers, however, due to its insufficient description of fine-scale scattering properties of needle shoots representing the main foliage elements. Therefore, we hypothesize that the second manner of constructing tree representations, in this case from user-defined 3D objects (Figure 1b), might be more appropriate for conifers and allow for producing more accurate DART simulations. The approach imports trees as groups of geometrically explicit 3D objects with different optical properties (i.e., trunks, branches, needle shoots of different age classes, etc.). DART can be executed directly over the original 3D objects constructed from triangular facets, which represents structurally the most accurate solution.
However, this option can be computationally highly demanding, especially when the imported 3D objects contains a lot of triangular facets. A third option is to transform the imported 3D foliage objects by voxelating them into leaf turbid cells (Figure 1c). Although such transformed crown representations are structurally less precise, they preserve all original aspects of crown shapes (e.g., distribution of air gaps) and foliage distributions (i.e., LAI and LAD). The voxelated nature of a forest stand representation substantially accelerates the computation and subsequently makes DART more operational with regards to its inverse modelling. Accurate physical representations that can account for complex 3D structures of ecosystems and link RS observations with studies of natural ecosystem processes are highly desirable. There is increasing demand, therefore, also for realistic 3D plant representations with natural architecture and foliage distribution. One group of methods for constructing 3D plants uses a set of parametric growing rules related to plant anatomy, topology, growth, and environmental factors derived from empirical observations. These are grammar-based methods, such as L-Systems [31], Grogra [32], OpenAlea [33], and recently developing computer graphic software (e.g., OnyxTree, Xfrog, and many others summarized at http://vterrain.org/Plants/plantsw.html). For instance, Widlowski et al. [34] used different virtual forests generated with Xforg and Arbaro software for the fourth phase of the radiative transfer model intercomparison (RAMI-IV) exercise. Nevertheless, most software tools create trees suitable only for exterior visualization and simulation of virtual landscapes (e.g., for computer games). They are less suitable for ecosystem studies requiring "biologically" correct 3D plant representations that account also for irregularities caused by natural and anthropogenic factors (e.g., spruce needle defoliation and regeneration; [35]).
Another group of methods for constructing 3D tree representations uses the rapidly developing technology of terrestrial laser scanning (TLS, also known as terrestrial LiDAR), which captures well the explicit 3D vegetation architecture [36][37][38]. Several studies have reported reconstruction of tree wooden skeleton (stem and branches) using TLS data by constructing a quantitative structure model (QSM) of topologically and hierarchically ordered components [39][40][41]. The TLS-based wooden skeleton reconstruction prefers data collected during leaf-off conditions, which means it is less than perfectly suited to conifers and evergreen-broadleaf mixed forests. In cases of evergreens, a TLS point cloud must be separated into wooden and foliage parts using (i) geometric properties [42], (ii) radiometric properties [43] or (iii) a combination of the two [44].
Considerably fewer studies have reconstructed a whole tree by adding foliage elements to a wooden skeleton. These studies concern mainly broadleaf tree species [45][46][47], and very few such attempts have been made for conifers [48]. Statistical information about spatial density, size distribution, and angular orientation of flat leaves could be extracted from the point cloud data [46,49,50]. For instance, Åkerblom et al. [46] introduced an algorithm to insert primarily broadleaf geometric 3D representations into QSM, but its extension to complex needle leaf objects considerably increased the computation requirements. Côté et al. [48] introduced the L-Architect model for conifers to reconstruct branching structure and to add foliage based on light availability derived from the point cloud data. Nonetheless, addressing 3D structure of conifers remains a challenge due to difficulties involving (i) leaf-wood separation, (ii) self-occlusion effects, and (iii) the complex geometrical nature of needle shoot elements.
The objectives of this study, therefore, were to (i) design an algorithm for reconstructing a close-to-reality 3D representation of Norway spruce trees from TLS measurements (i.e., reconstruction of wooden skeleton and natural distribution of needle foliage) and (ii) evaluate the effect of canopy structure on forest reflectance as simulated in a radiative transfer model. The second objective was achieved through DART simulations of four scenarios using different levels of spruce canopy structural details. The results of the scenarios were cross-compared and assessed against corresponding airborne hyperspectral and Sentinel-2 MSI RS reflectance observations.

Study Sites and Input Data
Computer reconstruction of trees, parameterization of DART, and the final assessment of the simulated canopy reflectance using airborne and satellite RS images required multiple data sources. All required data for this study were collected at two Norway spruce (Picea abies, L. Karst) forest sites located in the Czech Republic. The primary study site is located at the experimental research station Bílý Kříž in the Moravian-Silesian Beskydy Mountains (49°30 N, 18°32 E, 850-900 m a.s.l.). The site facilitates long-term research of tree ecophysiology and carbon fluxes. It is part of the national Czech Carbon Observation System (CzeCOS) and the international Integrated Carbon Observation System (ICOS) research networks e.g., [51,52]. The targeted stand is a spruce plantation growing on a southerly oriented slope (mean slope of 12.5°). Its characteristics in 2016 (when field and RS data were acquired) were as follow: mean tree age of 40 years, mean tree height of 17 m, mean stem DBH of 20 cm, mean stand density of 1230 trees ha −1 , and mean leaf area index (LAI) of 8.6 m 2 m −2 . The following additional stand characteristics were collected during a field and remote sensing campaign in August 2016 [53]: (i) foliar biochemistry and optical properties of needle samples for parameterization of DART simulations, and (ii) RS data (airborne hyperspectral and Sentinel-2 MSI images) used for evaluation of the DART reflectance simulations. Moreover, forest structural data on angular distribution of spruce shoots were obtained from an earlier study of Barták [54,55], who investigated 35-year-old spruce trees growing in the vicinity of the Bílý Kříž site.
Terrestrial laser scans (TLS) for 3D tree reconstruction could not be obtained for the Bílý Kříž site as the high density of forest stands did not allow collecting good TLS data on individual trees. Therefore, TLS data were acquired at a secondary study site, known asČerná hora and located in the Czech Republic's Šumava National Park (48°58 22.4 N, 13°32 53.8 E, 1250-1300 m a.s.l.). The spruce stand atČerná hora (100-year-old, mean height of 25 m, stand density of about 800 trees ha −1 , SSE-oriented slope with gradient of 7°) had been severely damaged by multiple bark beetle outbreaks just 1-2 years before the TLS data acquisition. Therefore it was possible to acquire high-density laser scans of surviving individual trees that previously had been growing inside a canopy.

Forest Structure Data
Three-dimensional reconstruction of virtual spruce trees required two types of forest structural data: (i) TLS data for reconstruction of woody skeleton and spatial distribution of foliage (clumping) and (ii) data on geometrical and angular distribution of spruce needle shoots.
The TLS data were acquired using an OPTECH Ilris-36D (Teledyne Optech, Vaughan, ON, Canada) scanner at theČerná hora site in October 2009. The Ilris-36D TLS emits laser beams at the wavelength of 1.5 µm with a maximum field-of-view window of 40°× 40°. It records the amount and intensity of returned laser pulses at up to 2500 points per second. The laser scanning data were acquired for 15 mature spruce trees that had been left standing after clearing trees damaged by multiple bark beetle outbreaks. Each tree was scanned from two opposed positions, recording the first and the last echo return with the point sampling distance of 2 cm. The point clouds of the first and last return from the two scanning positions were aligned and co-registered in a single geometric coordinate system using the PolyWorks IMAlign module (InnovMetric Software Inc., Quebec City, QC, Canada). The PolyWorks IMSurvey module was used to filter out individual tree point clouds and to separate wooden structures from foliage points based on the intensity of the reflected photon signals. Considering that green foliage reflects about 10% and woody elements up to 50% of the laser signal emitted at 1.5 µm, we determined manually the intensity thresholds between foliage and wooden elements for each scan (first and last returns from two opposite scanning positions) separately. The intensity thresholds varied between 8% and 18% depending on a scan. The extracted point cloud of wooden elements was subsequently used to reconstruct a 3D tree skeleton consisting of trunk and main branches, whereas the foliage point cloud was used in the algorithm for biologically correct distribution of needle shoots within a tree crown (as described in Section 2.2).
It was too difficult to derive the angular distribution of coniferous shoots within a crown from the TLS data. Therefore, we extracted data from a destructive study of Barták [54,55], who had measured the inclination of spruce shoots in respect to their age and position within a crown (see Figure 2). Field sampling took place at the Bílý Kříž site in 1990-1991 and twenty 35-year-old spruce trees were analysed and we assumed that the measured inclination angles are representative for the current forest stand of the similar age.

Needle Optical Properties
Optical properties (i.e., directional-hemispherical reflectance and transmittance) of spruce needles were measured during a field campaign at the Bílý Kříž site in August 2016. Needle samples of different age categories (current-, 1-and 2-year old needles) were collected from a sunlit and shaded branch of three trees growing at the main study plot in the vicinity of the flux tower. Needle optical properties were measured using a FieldSpec 4 spectroradiometer (Analytical Spectral Devices, Boulder, CO, USA) combined with an ASD RTS-3ZC integrating sphere using the method of gap fraction corrections suggested by Mesarch et al. [56] and further adjusted by Malenovský et al. [57] and Yanez-Rausell et al. [58]. Average needle spectra were used to parameterize the forest virtual scenes in DART, as described in Section 2.3.

Remote Sensing Data
Two types of RS data-airborne hyperspectral images and Sentinel-2 satellite data-were used for verification of DART-simulated canopy reflectance. Both RS images were acquired over the Bílý Kříž study site on 31 August 2016, close to the solar noon. Airborne hyperspectral images were acquired with a CASI sensor (Itres Research, Calgary, AB, Canada) in the range of 372-1044 nm (72 bands with spectral sampling distance of 10 nm and pixel size of 1.0 m) and with a SASI sensor (Itres, Research, Calgary, AB, Canada) in the range of 957-2442 nm (100 bands with spectral sampling distance of 15 nm and pixel size of 2.5 m). Corrections of the hyperspectral images were carried out in a data processing chain established at the CzechGlobe institute [59]. Radiometric corrections were performed using the factory calibration coefficients in the post-processing software developed by the sensors' producer. Geometric corrections (i.e., image orthorectification and geo-referencing) were performed using the GeoCorr software provided also by the sensors' manufacturer. Atmospheric and topographic corrections were made using ATCOR-4 software [60]. Images from CASI and SASI were combined into a single image hypercube (372-2442 nm) with final pixel size of 2.5 m.
The Sentinel-2 MSI image was corrected for atmospheric and topographic effects using the Sen2Corr toolbox provided by the European Space Agency. Spectral bands with nominal pixel size of 60 m and wide spectral bands were discarded, so that only the B2 to B7, B8a, B11, and B12 bands were retained in the final multispectral image cube.
Mean canopy hemispherical-conical reflectance function and its standard deviation were extracted from the CASI-SASI and Sentinel-2 MSI images for a region of interest (100 × 100 m) over the spruce forest study site in the vicinity of the flux tower, where field samples for needle optical properties and LAI measurements were obtained. The size of the region of interest was chosen to cover at minimum 4 × 4 Sentinel-2 pure forest pixels.

Reconstruction of 3D Spruce Tree Representation
The TLS point clouds supplied the explicit information about the position and orientation of dominant woody structures (trunk and branches) as well as about the foliage distribution and density of individual scanned trees. Woody and foliage elements were separated first using an intensity threshold that was determined for each tree scan manually based on visualization of laser return intensity histograms. The TLS data was itself insufficient for deriving the angular distribution of spruce needle shoots (considered here as the elementary units). Therefore, this information was taken from an external study on Norway spruce tree architecture [54,55] Figure 2.
Reconstruction of a 3D spruce representation from the TLS point cloud data was executed in the three stages documented in Figure 3:

1.
Spatial separation and reconstruction of the wooden skeleton including trunk and main branches using an automated algorithm developed by Sloup [41] (Section 2.2.1, Figure 3b,c).

2.
Translation and scaling of the foliage point cloud and the reconstructed wooden skeleton to fit site-specific tree dimensions. As the TLS data had been acquired at theČerná hora site for taller trees, a scaling of dimensions was required to fit tree sizes at our Bílý Kříž primary research site (Section 2.2.2, Figure 3d).

3.
Distribution of shoots within the tree crown by combining the information about foliage density from the TLS data with the ancillary information on shoot angular distribution [54,55]. (Section 2.2.3, Figure 3e).

Reconstruction of Wooden Skeleton-Trunk and Main Branches
The algorithm for reconstruction of tree wooden skeleton had been developed for Norway spruce trees by Sloup [41] using TLS data acquired at theČerná hora site. It uses the method of Verroust and Lazarus [61] for extracting skeletal curves. It is fully automated and able to process point cloud data containing gaps from the omission of laser returns due to shadowing of any obstacles. It identifies spatially-related clusters of points (components) and the process of connecting branch components preserves the real architecture of branches and their connections to the trunk. The algorithm reconstructs the trunk and main branches, small branches at the crown periphery could not be reconstructed, because of insufficient details captured by TLS scans. Small wooden parts such as twigs can be included in the shoot representation.
The algorithm can be summarized into the following three steps (details and the code can be obtained from https://is.muni.cz/th/325196/fi_m/?lang=en):

1.
Component identification: spatially related clusters of TLS points are identified and separated.
constructing a neighbourhood graph to connect close points [61]; the result is a set of components, each consisting of spatially related points belonging to the same branch or trunk part.

2.
Component analysis: branch-like structures are reconstructed for each identified component.
constructing a geodesic graph for each component that describes how the distance from the designated source point gradually increases throughout the component [61]; -collapsing points with similar distance from the source point, unless they belong to a different sub-branch.

3.
Component connection: components are interconnected to form final trunk-branch structures.
initial estimation of a trunk by fitting a segmented conical model to automatically selected trunk components; -iterative selection of components to be connected to any other component while minimizing a cost function (combining the euclidean distance of the endpoints, outgoing vectors describing directions and probabilities of branch continuations, and slope function that roughly estimated expected slope of the main branches).

Translation and Scaling of Foliage Point Cloud
A reconstructed 3D wooden skeleton was built in a local coordinate system with its origin at the centre of the trunk base and neglecting original coordinate system of point cloud. A translation matrix between the two was established to ensure the same coordinate system for the foliage TLS point cloud as for the wooden skeleton. The reconstructed wooden skeleton and foliage point cloud were subsequently scaled to match site-specific tree dimensions. For this purpose, a scaling factor, f s , was defined as the ratio between the required tree height and the maximum z-coordinate of the reconstructed wooden skeleton.
In our case study, the TLS data had been obtained from trees older than the spruces at our primary research site. Despite the difference in tree age of almost 60 years between theČerná hora and Bílý Kříž study sites, the discrepancy in the mean canopy height is not so very pronounced, as both forests are growing in mountainous environments and are exposed to similar environmental stress impacts. We therefore assume their woody and foliage architectures to be similar and possibly transferrable from one site to the another one. Scaled tree dimensions (crown diameter and height) were compared with existing forest allometry data showing good agreement. In addition, only foliage distribution is derived from the TLS point cloud data, not the actual density and angular distribution (described in the next section).

Shoot Distribution within a Crown
The third stage of tree reconstruction consists in distribution of the main foliage elements (i.e., needle shoots) within a tree crown. Because biochemical compounds together with the spatial and angular distribution of shoots influence light absorption and scattering, we considered in our reconstruction two needle age categories: (i) current-year needles, growing mainly at the crown periphery, and (ii) older needles, growing mainly inside the crown. The shoot distributing algorithm uses information about shoot positions from the foliage point cloud (a higher density of laser points indicates a higher density of leaves). The actual shoot density is controlled by the user through the definition of tree LAI. The leaf angle distribution (LAD) was adopted from the measurements of Barták [54] Figure 2. That algorithm has three main steps: (i) calculation of shoot positions, (ii) division of shoots into two age classes, and (iii) placement of shoots at predefined positions and angular orientations.

Calculation of Shoot Position
The number of shoots is given by the user-defined tree LAI. The total number of shoots within a reconstructed tree crown is computed as: where n shoots is the total number of shoots within a crown, S shoot is total area of needles in one shoot (e.g., as calculated in the Blender software, www.blender.org), and S proj is the circular projection area of the tree crown represented by the extracted and scaled foliage point cloud. The diameter d of the circular projection area S proj is calculated as: where max x ,max y , min x , and min y are the maximum and minimum coordinates of the foliage points in the x and y dimensions. Calculation of the shoot positions uses a grid of voxels. The foliage point cloud is, therefore, first separated into a grid of voxels of a given size (in our case 0.18 m). The voxel size influences the computational time, as placing more points within a single voxel requires more time for computing all shoot positions. Hence, the voxel size should be user-specified according to the computational resources available and the density of the foliage point cloud.
The actual shoot positions are calculated using the k-means function based on the cluster analysis theory. Each voxel contains a set of points x = {x 1 , x 2 , . . . , x n }, where n ∈ N is the number of points inside the voxel currently being processed. Each point is defined by its position x i = (x, y, z). The k-mean clustering partitions points into k < n subsets S = {S 1 , S 2 , . . . , S k }. The k number of subsets is then computed from the total number of shoots within the crown (n shoots , Equation (1)) as: where n points is the total number of points in the whole foliage point cloud. The k-mean function computes also the coordinates of the cluster centroid, which then represents the actual shoot position. The algorithm iterates through all voxels containing at least one laser point.

Separation of Shoots by Age Categories
For this step, we need to identify a tree crown envelope to assign current shoots at the periphery and older shoots inside the crown. Therefore the foliage cloud point has to be divided into a voxel grid again and we recommend setting the size of voxels at approximately twice larger than that for the shoot positioning grid (in our case we use the voxel size of 0.3 m). The tree crown envelope and its width were defined using an iterative approach that is presented in Appendix A. The tree was then divided into several vertical layers (l n ) and a percentage of current-year shoots was assigned to each layer based on the field measurements, as summarized in Table 1 (details in Appendix B). Different optical properties ( Figure 4) were defined for each needle category in order to capture the natural biochemical variability.   Table 1. Field data on spatial and angular distribution of Norway spruce shoots (destructive sampling of 20 trees) extracted from Barták [54,55]. Percentage of current-year shoots, mean shoot elevation angle, and that angle's variation in the vertical crown levels associated with branch whorls are counted from top to base. Because needle shoots can be constructed as 3D objects, their shape and complexity can be altered according to user and computational requirements. Here, we used two different 3D shoot models: (i) a shoot constructed from individual non-flat needles along a central twig (Figure 5a), and (ii) a shoot simplified to radial planes with total area equal to the leaf area of single needles while maintaining the original needle lengths (Figure 5b). The main difference between the shoot models consists in the number of object facets. The shoot model with individual needles contained 824 facets per shoot, whereas the simplified model with planes contained only 16 object facets. Both shoot representations had been created and used for the RAMI-IV radiation transfer model intercomparison exercise [34]. Once the exact positions for shoots have been determined, shoots as individual 3D objects are distributed within the tree crown. The shoots must first be rotated according to predefined angular boundaries varying in the vertical direction as summarized in Table 1 and Figure 2. Finally, shoots are translated to their actual positions inside the crown. A transformation matrix for angular rotation and spatial translation is calculated for each shoot separately, as described in Appendix C.

DART Forest Scene Parametrization
Small forest scenes of 10 × 10 m and with a fixed distribution of 10 reconstructed trees were created in DART using either internal basic tree shapes filled with a turbid medium (Figure 1a, Appendix D) or the external 3D representations reconstructed from TLS point clouds (Figure 1b,c). Main differences among the parameterization approaches are not only in the tree structure descriptions but also in definitions of the leaf optical properties, LAD, and LAI. The remaining DART input parameters were tailored to forest conditions captured at the Bílý Kříž study site during the field and remote sensing campaign in 2016. These include tree dimensions; optical properties of forest floor, spruce bark and needles; actual sun position during the acquisition of airborne and Sentinel-2 images; and description of the atmosphere (based a on mid-latitude summer model in Modtran [62]). The input parameters common to all DART simulations are summarized in Table 2, and the optical properties are depicted in Figure 4c. eaf optical properties were computed differently for the trees created internally and externally. In the case of basic DART tree shapes, they were mixed based on the ratio of current-year and older shoots (Figure 4a), as follows: (i) 7:3 for a sunlit vertical crown layer (top 17% of tree crown), (ii) 1:1 for a transitional sunlit and shaded crown layer (middle 30% of tree crown), and (iii) 1:8 for a shaded crown layer (bottom 53% of tree crown). In the case of the reconstructed 3D trees, the leaf optical properties were assigned according to the distribution of two needle age categories as described in Section 2.2.3 (Figure 4b).

Evaluation of the Spruce Forest RT Simulations
All RTM simulations were conducted with DART model version 5.7.1 (build 1061) according to four scenarios representing various levels of canopy structural complexity as summarized in Table 3. Scenario A, basic turbid, which is regarded as structurally the simplest, uses basic predefined crown shapes and represents the foliage as turbid cells. The remaining three scenarios, B-D, use 3D trees reconstructed from TLS, but they differ in the complexity of 3D shoot objects and in the foliage representation as object facets or as turbid cells transformed from these objects. Scenarios A-C were simulated for the spectral range from 450 to 2300 nm with the 2 nm bandwidth and subsequently convolved to the spectral resolutions of CASI-SASI and Sentinel-2 MSI images. Scenario D, needle facet, is structurally the most complex and thus computationally the most demanding scenario. Therefore, only 13 spectral bands with the central wavelengths of Sentinel-2 MSI bands were computed for this scenario. From DART-simulated images of each scenario we extracted mean scene spectral reflectance and its standard deviation, as well as mean reflectance of sunlit crown parts and its standard deviation. DART-simulated spectral signatures were compared with the canopy reflectance signatures extracted from CASI-SASI and Sentinel-2 images. Root mean square error (RMSE) and index of agreement d between simulated (R SI M ) and observed (R OBS ) spectral signatures were computed as: where n is total number of either CASI-SASI or Sentinel-2 MSI spectral bands, R OBS is mean value of observed spectral signatures. The CASI-SASI spectral bands that are strongly influenced by atmospheric water vapour absorption (i.e., 893-1002 nm, 1107-1167 nm, 1302-1527 nm, and 1737-2052 nm) were removed from the RMSE and d computation. Further, we extracted and cross-compared multi-directional reflectance distributions of the simulated canopy scenarios.
The main criteria to evaluate the optimum scenario are (i) overall agreement with real observations extracted from airborne hyperspectral and Sentinel-2 image data and (ii) computation requirements.

Reconstructed 3D Spruce Trees
The total number of leaf facets of a single tree constructed according to TLS scans ( Figure 6) varied between 11 million and 30 million, depending on crown dimensions (see Table 4). The LAI of each individual tree crown, i.e., relative to the tree crown vertical projection onto the ground, was set to 12 m 2 m −2 (Equation (1)). This value had been found to be the optimal LAI for building up a forest scene with canopy cover of approximately 80% and canopy LAI equal to 8.5 m 2 m −2 (as measured in the field). The DART forest scene of 10 trees was built from the four individual tree representations, which were replicated, rotated, and re-distributed to achieve the desired canopy cover while avoiding a regular spacing (Figure 7). Table 4. Structural characteristics of the four reconstructed 3D spruce trees ( Figure 6) populated with complex needle shoot objects (leaf area index = 12 m 2 m −2 ). Rigorous validation of the reconstructed 3D trees would be difficult and tedious, as it would require that the evaluated structural parameters of the reference trees be measured independently from their TLS scans. Moreover, the original TLS data were spatially scaled to different dimensions and the density of needle shoots was altered according to the user-defined value of tree crown LAI. Therefore, the distribution of shoots and the overall habitus of reconstructed trees were inspected only through visual comparison with the original TLS point clouds (see Figure 3a) and photographs of modelled trees. Later, we evaluated the reconstruction of 3D trees in the context of the forest stand canopy reflectance simulated in DART and compared to the real airborne and space-borne observations.

Simulations of Forest Canopy Reflectance in DART
DART simulations of Norway spruce forest canopy reflectance between 450 and 2300 nm (i.e., 926 bands with step and width of 2 nm) for scenarios A-C were carried out on a local server (40 cores, 350 GB RAM). Simulations of the turbid cell scenarios A and B were computed relatively quickly (each scene in less than 3 h, using only 2 processor cores), but to retain the foliage elements as object facets required substantially increased computation time. The forest scene of scenario C (shoots as planes) contained 3.7 million leaf facets and its computation took almost 34 h (using 8 processor cores). Computation of scenario D, needle shoots, containing 189 million facets, was extremely demanding computationally and had to be processed on the Czech national computer grid MetaCentrum (https: //metavo.metacentrum.cz/en/index.html). Figure 7 displays the reconstructed virtual forest scenes and their corresponding false-colour composite images of canopy reflectance simulated in DART. Scenario A, basic turbid, shows regular shapes and homogeneous filling of spruce crowns, which are not often seen on real airborne spectral images. The structurally more heterogeneous scenarios B-D revealed larger spatial variability of reflectance within and outside of the crown pixels. One can observe more shaded pixels inside and at the crown edges and also more sunlit forest floor pixels (the cyan colour), especially in scenarios C and D. Despite that crown diameters were the same in all scenarios, retaining the real branch architecture resulted in tree crown peripheries being more transparent when compared to scenario A.
Overall, the greatest simulated canopy reflectance in green and near-infrared (NIR) spectral domains was found in scenario A. Simulated red and shortwave infrared (SWIR) reflectance values were converging in all simulated scenarios, which indicates that these spectral wavelengths are less affected by differences in simulated canopy architectures (Figures 8). All simulations were carried out with added sensor viewing directions around the hotspot (i.e., hotspot oversampling). One can observe that turbid media (scenarios A and B) produced a sharper hotspot feature compared to facet-based simulations (scenarios C and D) (Figure 8), which is an artifact of the turbid canopy hotspot analytical computation. Inspecting closely the multi-angular reflectance behavior of the four scenarios (Figure 9), one can note overestimations in the green wavelength (552 nm) for the turbid scenarios (A and B) when compared to the facet-based scenarios (C and D), especially at oblique zenith viewing angles. All multi-angular behaviors agree well for the red reflectance band (666 nm), as they were influenced by strong chlorophyll absorption. The only exception relates to the plane facet of scenario C, which indicates a reflectance overestimation in oblique back-scattering directions. The largest multi-angular reflectance differences can be seen in the NIR band (776 nm). The basic turbid scenario A overestimates, while the plane facet scenario C underestimates, NIR reflectance in all viewing directions when compared to scenario D. Finally, all simulated scenarios produced a consistent reflectance with minimal differences in the SWIR wavelength (1647 nm). Minor overestimations were observed for the needle turbid scenario B in strongly oblique (> 60°) forward-scattering viewing directions.
The first derivative (FD) of simulated canopy reflectance reveals changes in the shape of spectral signatures. FD was not calculated for scenario D, for which shape changes would be too abrupt, due to its inclusion of just 13 simulated spectral bands. Results for scenarios A-C show that, despite their differences in absolute reflectance, the signature shapes are very similar, especially for the NIR and SWIR wavelengths that are not impacted by atmospheric water vapour absorption ( Figure 10). The largest differences were found in the green and mainly in the red-edge region. One can observe substantial FD peaks, which correspond to inflection points on the steepest parts of the reflectance curves. In most bands, the FD of scenarios A, basic turbid, and B, needle turbid, are similar, which likely is caused by their turbid cell nature. Interestingly, the reflectance signature shape closest to that of the real airborne data acquired over the study site using CASI and SASI sensors was found in the case of scenario C, plane facet.

Comparison with Airborne and Satellite Data
Simulated canopy reflectance signatures for the four scenarios (A-D) were plotted and compared against spectral signatures extracted from the corresponding airborne hyperspectral CASI-SASI images and satellite Sentinel-2 MSI image (Figures 11 and 12). Taking advantage of the high ground spatial detail of airborne hyperspectral images, we compared the simulated scenarios with (i) mean forest stand spectra extracted from CASI-SASI image with pixel size of 2.5 m and entire DART scenes ( Figure 11a,b), (ii) spectral signatures extracted from sunlit tree crowns of CASI and DART images of 1.0 m pixel size (without contributions from canopy shaded parts and background) (Figure 11c,d), and (iii) mean canopy spectral signatures extracted from Sentinel-2 MSI image with pixel size of 20 m and DART scenes (Figure 11e,f). The comparison of forest stand reflectance extracted from airborne CASI-SASI and Sentinel-2 MSI images for the same region of interest showed good agreement, even though the reflectance variability in the airborne data is naturally greater due to its finer pixel size (2.5 m) (cf., Figure 11a,e). The consistency between the two independent RS observations gives us a certain level of confidence in using the RS data as a robust benchmark for evaluating the different scenarios of spruce forest parameterization in DART. When compared to the RS observations, all scenarios (A-D) overestimated the stand canopy reflectance across the entire spectral region between 450 and 2300 nm. Scenario A, basic turbid, yielded the largest differences from the RS observations, mainly in the visible and NIR regions. Other implementations of the tree 3D architectures (scenarios B-D) brought the DART forest reflectance simulations closer to the RS observations.
In the visible region, all scenarios systematically yielded greater reflectance than did RS observations. The simulations manifest a systematic reflectance overestimation by 0.01-0.02 in green peak. The red reflectance of the chlorophyll-a absorption feature around 675 nm also exceeded that of the RS observation by 0.01, and in the case of scenario C, plane facet, this was by even 0.02. This deviation disappeared, however, when inspecting the reflectance of only sunlit parts of tree crowns (cf., Figure 11b,d). This suggests that the mean reflectance of scenario C might be influenced by sunlit background pixels that are substantially increasing simulated values in red wavelengths. The largest reflectance variability among the simulated scenarios and differences from RS data were obtained for the NIR spectral region. When compared to the airborne data, the NIR reflectance around 850 nm of scenario A, simple turbid, was as much as 0.1 larger. This difference increased to 0.15 for sunlit crown spectra. A similar NIR reflectance overestimation was observed also for scenario D, needle facet. NIR reflectance of scenarios B, needle turbid, and C, plane facet, fell within the standard deviation interval of the RS observations. Finally, all scenarios simulated very similar reflectance for the SWIR region. Nevertheless, all scenarios systematically overestimated the RS observations around 1650 nm by 0.04 and around 2200 nm by 0.02. Figure 12. Scatterplots between DART-simulated and remote sensing (RS) data extracted from spruce forest stand reflectance. Top row: average stand reflectance extracted from CASI-SASI image, middle row: average reflectance of sunlit tree crowns extracted from CASI image, and bottom row: canopy reflectance extracted from Sentinel-2 MSI image (only narrow spectral bands with spatial resolution ≤20 m that are not affected by atmospheric water absorption are plotted). In order to provide statistically comparable data sets, the open circles in the case of airborne data (top and middle row) represent wavelengths corresponding to the Sentinel-2 MSI spectral bands.
Scatterplots, similarity indicator d, and RMSE between the DART-simulated scenarios and RS observations (airborne and Sentinel-2 MSI) were plotted and computed only for wavelengths corresponding to those Sentinel-2 MSI spectral bands not affected by atmospheric water absorption, (i.e., not interpolated during the atmospheric corrections) (Figure 12). The highest d and lowest RMSE values were obtained for scenario C, plane facet, followed by B, needle turbid, and D, needle facet.

Reconstructed 3D Spruce Trees
The first goal of this study was to create a biologically correct 3D digital representation of Norway spruce trees for purposes of forest canopy radiative transfer modelling [63,64]. The results demonstrate that the designed algorithm is capable of reconstructing the architecture and geometrical properties of spruce wooden and crown foliage elements from TLS observations of individual trees. Quantitative validation of the reconstructed 3D spruce representations could not be performed, as doing so would require reference measurements of tree structural (geometrical) features and the scanning device position, which were not at our disposal. Such a detailed quality assessment is not even needed, however, because the aim was not to reconstruct a given tree as accurately as possible but rather to obtain biologically correct 3D tree objects that can form a representative forest canopy with a variable leaf area index and user-defined foliage biochemical traits. Consequently, the developed method allows one for reconstructing trees with any required value of LAI in order to match the forest stand LAI measured at the site where the RS data was acquired.
The algorithm is automated and works with separated point clouds of wooden skeleton and foliage elements. The fidelity of a reconstructed tree architecture depends on the quality of TLS scans, which drives the decision on either manual or automatic separation of wood and foliage point clouds. The wooden skeleton reconstruction algorithm [41] correctly handles point gaps in wood point cloud data caused by occlusion of emitted laser beams due to branches within the scanner field of view. The algorithm for foliage distribution has been designed specifically for Norway spruce trees, but it might be applicable also to other spruce species and, after necessary adjustments, to other tree biological species. A user has the possibility to control tree LAI by adjusting the density of foliage objects (in this case shoots) in proportion to the density of TLS points. The angular orientation of spruce shoots had to be set according to manual field measurements of a limited sample of trees [54], as it was not feasible to retrieve this information from TLS scans. The leaf angle distribution might be obtained for broadleaf trees e.g., [49], but the method is not applicable for conifers with small narrow needles organized in shoots.
The approach developed in this study can realistically reconstruct the architecture of an entire spruce tree, including wooden skeleton and foliage geometry, while also implementing different groups of needle shoots according to their age. This ability is crucial for such coniferous tree species as Norway spruce, where leaves or shoots emerging in different vegetation seasons (i.e., several generations) co-exist within a single tree crown. Previous studies (e.g., [63][64][65][66]) solved this problem by dividing tree crowns into several vertical layers, to which different needle optical properties are attributed according to measured needle-age ratios and sun radiation availability. Our approach allows one for populating a tree skeleton with real needle shoots, thus reflecting the realistic distribution of specific leaf optical properties according to their age. This capability is especially important for simulations of high spatial resolution optical images (with pixel-size on the order of centimetres), where diversity in both, optical properties and forest structure has an important influence on the forest canopy reflectance. In this respect, the reconstruction of wooden skeleton and crown architecture might be sufficient, but the distribution of optical properties could be treated more realistically. Besides shoots, young branches and younger parts of the trunk are in the case of spruce also populated by needles. In our study, however, all branches and trunk parts were assigned with the optical properties of bark. Implementing needles, or at least their optical properties, for the top parts of branches and trunks might mean a substantial difference in simulated red and NIR canopy reflectance, especially for the nadir viewing direction.
Another possible improvement is to create a more realistic reference needle shoot 3D model. The representation used in this study was relatively simple. Its design was ruled by the limited capabilities of handling 3D objects in older versions of DART and limited availability of computation resources. Moreover, a universal distribution of needles along the central wooden twig does not resemble the actual variability existing in natural shoot geometrical formations. The architecture of each shoot varies depending on illumination conditions, age, and stress impacts. More genuine 3D shoot representations could be obtained from high-resolution laboratory laser scanning. Implementing several types of shoot representations reflecting the natural variability driven by shoot age and location within the tree would improve the accuracy of the radiative transfer modelling, and particularly that of the light scattering within shoots.

Simulations of Forest Canopy Reflectance in DART
One of the key issues for remote sensing of coniferous forests is the upscaling of optical properties from needles to canopies. From an operational point of view, it is easier to measure optical properties of needles rather than of shoots. Inasmuch as shoots are the main scattering elements in coniferous canopies, however, there is a strong demand for characterizing the shoot optical properties and phase functions with high accuracy [67][68][69][70]. Traditional forest reflectance models based on the radiative transfer equations handle the shoot-level needle clumping by correcting the radiation attenuation coefficient with a clumping index. That compensates for a reduction in the interception of radiation by the canopy at a fixed LAI [68]. Hence, implementation of the precise 3D spruce representations can help in disentangling and better understanding light interaction related to the needle foliage clumping at different structural levels (i.e., shoots/twigs, branches, crowns, and stands).
The reflectance simulation of whole forest canopy in DART served to compare different approaches of tree parametrizations and to assess their performance in terms of accuracy, operational implementation, and feasibility. To produce DART simulation of a forest canopy scene 10 × 10 m in size and with 3D tree objects created from almost 189 million facets was until recently unfeasible. We were able to execute such a geometrically detailed simulation using high-performance computation resources of the Czech Republic's national grid organization known as MetaCentrum. Even then, we were able to simulate only 13 spectral bands and total computation time took more than three weeks. Outputs of this most-detailed simulation were used in comparison with other, simpler parametrizations, where we are able to simulate hyperspectral images with the 2 nm spectral resolution in much shorter time and using substantially fewer computational resources. With future DART code improvements and availability of more powerful computation resources, we expect to simulate full hyperspectral images even for the most detailed 3D representations, as in the scenario D.
The results confirmed an expectation that simulations using the reconstructed 3D representations produce more consistent canopy reflectance, especially in red and NIR wavelengths (see Figure 8). The first derivative of the canopy reflectance revealed that the most substantial shape changes in the reflectance function are between the basic turbid scenario A and the other scenarios at the red edge spectral region ( Figure 10). Nevertheless, depending in each case upon the purpose of their use, the results of all scenarios could be found acceptable. Scenario A, basic turbid, differs most from the others. It is the least demanding of computational resources and its reflectance results are in agreement with the reference scenario D in the red and SWIR wavelengths. Contrary to this scenario D, needle facet requires high-performance computational resources and a substantial simulation time, due to detailed 3D needle shoots. Reflectance results of scenario B needle turbid are in most cases very similar to the most detailed scenario D needle facet. The only substantial difference was found in the NIR region (800-1200 nm). As expected, canopy reflectance of both scenarios simulating canopy as a turbid medium behaved differently in the hot spot region, due to the turbid simulation's hot-spot analytical solution. Finally, compared to scenario D, the scenario C plane facet slightly underestimated red reflectance, especially around the hot-spot region, and overestimated NIR reflectance in all simulated viewing directions (see Figure 9).

Comparison with Airborne and Satellite Data
Widlowski et al. [34] pointed out that the reconstruction of a natural forest stand scene for purposes of modelling radiative transfer could be problematic, because several tree habitus characteristics might be unavailable or difficult to measure. Many canopy biochemical, optical, and mainly structural properties, such as the geometrical distribution of foliage elements (leaves or shoots), in a canopy vertical profile or allometric parameters of wooden elements, including the branching angles, are often not measured simultaneously with RS spectral image acquisitions. The Bílý Kříž research station in Czech Republic is an ecological monitoring site where ground measurements of the monitored spruce stand were collected simultaneously with airborne imaging spectroscopy data (CASI and SASI sensors) and satellite Sentinel-2 MSI images. This allowed us not only precisely to parameterize the DART 3D scene, but also to compare the simulations with RS data acquired during the same vegetation phenophase and at the same location.
The comparison revealed a systematic overestimation of simulated bidirectional reflectance factors, especially within the spectral interval of 450-570 nm. This mismatch could be caused by a greater contribution of background, represented in our simulations by needle litter having relatively high visible spectral region of reflectance. The possibility that this is the cause is supported by results obtained from separating spectral responses of sunlit crowns within the simulated and acquired images. Suppression of the understorey influence in red wavelengths is obvious in the case of scenario C, plane facet, (cf. Figure 11b,d). Another possible cause of this discrepancy might be a mismatch in the distribution and content of shadowed pixels (between tree crowns) in RS data and simulations. The physical modelling treats light transfer according to optical laws without any noise, while the RS sensors and processing of their data introduce several types of noise that DART does not simulate. Although we used the needle optical properties measured during the RS data acquisitions, a certain level of non-representativeness originating from a disproportionate ratio between current and older needle generations might still occur and contribute also to this mismatch. Finally, a spectral shift between simulated and RS data around 1800 nm is caused by the interpolation of airborne data reflectance in the course of the atmospheric corrections.
The comparison of canopy reflectance extracted from RS data and DART simulations, as if mimicking narrow spectral bands of Sentinel-2 MSI with the spatial resolution ≤20 m and unaffected by atmospheric water vapour absorption (Figure 12), revealed almost the same results for the last three scenarios (B, C, and D). This suggests that the most efficient and optimal approach for simulating the Sentinel-2 MSI forest stand reflectance in DART is to use scenario B, needle turbid, which can deliver comparable reflectance simulations as scenarios C and D, but the computation time is substantially shorter.

Conclusions
In this study, we evaluated the influence of Norway spruce forest structural representations on DART modelled top-of-canopy reflectance. DART simulated nadir reflectance values were compared with actual airborne hyperspectral and satellite multispectral images acquired over the simulated research site. Because the multi-angular directional reflectance measurements of the investigated spruce forest were not available, the multi-angular simulations were compared against the most detailed and presumably most accurate scenario D, needle facet. The results revealed improvements in reflectance accuracy simulated for scenarios wherein the 3D forest scene representations were retrieved and parametrized using the method designed in this study. We identified scenario B, needle turbid, as the most effective scenario for nadir simulations. This scenario was able to produce the visible and NIR reflectance with an acceptable level of accuracy (i.e., similar to that of scenario D and close to actual observations) with relatively short computational time and without a need for high-performance computer resources. Being based on the needle turbid representation of the stand, however, scenario B could not very well reproduce reflectance around the hot-spot and in the extreme zenith angles. In the case of the multi-angular directional reflectance, scenario C, plane facet, tracks closely the most detailed scenario D, needle facet, while maintaining a comparable (in NIR) or even better match (in visible wavelengths) with the reflectance of RS data acquired from the nadir.
The 3D representation of Norway spruce trees reconstructed from TLS data with the semi-automatic algorithm designed in this study was found to reassemble adequately the habitus of real Norway spruce trees and increase the fidelity of canopy reflectance simulations in DART. The topographically accurate capturing of the tree structures and more detailed parametrization of foliage optical properties throughout a tree vertical profile allowed for more accurate solar radiation scattering within a simulated forest stand. Because the results obtained showed reflectance overestimations in visible and NIR spectral regions when compared with RS observations, there clearly is room for further improvements. Subsequent research should be focused upon (i) more realistic distribution of optical properties for wooden parts of forest stands (especially for young branches and tips of tree crowns), and (ii) the extension and improvement of needle shoot representations. Several needle shoot types could be reconstructed and used instead of a single, universal 3D model. Their size and geometrical parameters should reflect natural variability driven by shoot location within the tree crown, age, needle angular distribution, and within-shoot clumping in order to produce more realistic spruce tree representations for purposes of estimating radiative transfer.

Acknowledgments:
The authors would like to thank Petr Sloup and Tomáš Rebok for providing code for 3D tree wooden skeleton reconstruction and Jan Novotný for advice regarding terrestrial laser scanning data processing.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: iterate through every voxel containing at least one shoot position Figure A1. The process of finding crown envelope. The corresponding steps are described above 1.-4.
The crown envelope convexity reduction prevents new shoots from being placed inside the crown as they naturally grow at the periphery. This may happen when the air gap size is larger than the voxel size. The following steps were implemented to reduce the inappropriate crown envelope convexity (illustrated in Figure A2   iterate through every horizontal layer Figure A2. The process of crown envelope convexity corrections. This process describes in detail Step 4, consisting in finding the crown envelope. Corresponding steps 4.1. to 4.13 are presented above. Figure A3. The process of separating and finding crown envelope. The process iterates through every voxel, one by one, so long as there are still points from foliage point cloud. Figure A4. Extension of a crown edge towards the tree trunk. The extension is conducted in two consecutive steps, namely in the x and the y directions. Light green represents voxels contained in tree crown envelope, and dark green represents inside of the crown.

Appendix C. Angular Distribution of Shoots
The position in which the shoot was supposed to be translated is noted as L = (l x , l y , l z ). A shoot was first rotated around the y and z axis. The elevation angles β, listed in Table 1, are dependent on height (l z ). The calculation of the azimuth angle γ is based on l x and l y . It is randomized by an angle ψ getting values within the range − π 4 , π 4 (Equation (A1)).
The vector for shoot translation was calculated as v = (L − 0), where 0 is the origin of the coordinate system. The transformation matrix was computed as: Equation (A2) represents the general matrix for shoot models transformation, which was applied to distribute shoots of the two age categories within the reconstructed 3D model of a Norway spruce tree.

Appendix D. Parametrization of the Basic Turbid DART Tree Representation
Predefined tree shapes available in DART (Fugure 1a) allow constructing spruce trees as composed ellipsoids with precisely defined tree crown dimensions (Table A1). Each tree crown was further divided into three vertical levels. Within each level, different leaf optical properties were defined, tree stem diameter adjusted, and percentage of empty leaf cells (gaps) controlled. Gaps were distributed such that there were more gaps along the stem and at the bottom (which are shaded parts) and thus to mimic natural tree defoliation (Table A2). Leaf optical properties were computed as weighted means of leaf spectra obtained for needles of different categories (sun-exposed needles from the tree top vs. shade-adapted needles from the tree bottom and needles of different age categories ranging from the current year to needles up to 2 years old). Weighting factors (Table A3) were derived from data on distribution of needle age classes within spruce crowns obtained by field measurement conducted at the Bílý Kříž site by Pokorný and Marek [71]. Average leaf optical properties for each crown level are shown in Figure 4a. Leaf angle distribution was in this case fixed and defined as a spherical distribution. Forest scene leaf area index was fully controlled by DART and was set up to be equal to 8.5 m 2 m −2 . Table A1. Basic tree parameters used to construct spruce trees in DART using predefined tree shapes. The figure on the top illustrates the individual tree parameters (adapted from DART manual [30]).  Table A2. Definition of tree crown levels with inner defoliation zone that is illustrated by parameters a and b (the horizontal axis illustrates tree crown diameter, where 0 represents the tree stem and 1 represents the crown's periphery, adapted from DART manual [30]).