Simulating imaging spectroscopy in Tropical Forest with 3D Radiative Transfer Modeling

: Optical remote sensing can contribute to biodiversity monitoring and species composition mapping in tropical forests. Inferring ecological information from canopy reflectance is complex and data availability suitable to such a task is limiting, which makes simulation tools particularly important in this context. We explored the capability of the 3D radiative transfer model DART to simulate top of canopy reflectance acquired with airborne imaging spectroscopy in complex tropical forest, and to reproduce spectral dissimilarity within and among species, as well as species discrimination based on spectral information. We focused on two factors contributing to these canopy reflectance properties: the horizontal variability in leaf optical properties ( LOP ) and the fraction of non-photosynthetic vegetation ( NPVf ). The variability in LOP was induced by changes in leaf pigment content, and defined for each pixel based on a hybrid approach combining radiative transfer modeling and spectral indices. The influence of LOP variability on simulated reflectance was tested by considering variability at species, individual tree crown and pixel level. We incorporated NPVf into simulations following two approaches, either considering NPVf as a part of wood area density in each voxel or using leaf brown pigments. We validated the different scenarios by comparing simulated scenes with experimental airborne imaging spectroscopy using statistical metrics, spectral dissimilarity (within crowns, within species, and among species dissimilarity) and supervised classification for species discrimination. The simulation of NPVf based on leaf brown pigments resulted in the closest match between measured and simulated canopy reflectance. The definition of LOP at pixel level resulted in conservation of the spectral dissimilarity and expected performances for species discrimination. Our simulation framework could contribute to better understand performances for species discrimination and relationship between spectral variations and taxonomic and functional dimensions of biodiversity. DART model. We used ALS data in order to reconstruct 3D forest architecture.


Introduction
Tropical forest ecosystems host at least two-thirds of the terrestrial biodiversity and contain about 25 % of the carbon in the terrestrial biosphere [1,2]. The loss of biodiversity currently observed worldwide strongly impacts tropical forests [3]. This decline of biodiversity is induced by anthropogenic factors including land use change, exploitation of natural resources, and climate change, leading to forest disturbance and disappearance [4][5][6][7], with potential feedback due to climate-biosphere interactions [8].
Many studies have reported the relationship between ecosystem functions and biodiversity [9][10][11][12], others highlight ecosystem services provided by the biodiversity of such forest ecosystems and their importance sustaining life on Earth [13,14]. Thus in order to monitor biodiversity and ensure their conservation, a set of so-called Essential Biodiver-Numerous approaches have been implemented to model radiative transfer in terrestrial environments. One-dimensional models such as PROSAIL [43] have been widely used to simulate radiative transfer in spatially homogeneous canopies like crops, grasslands, and closed forests [44,45]. In complex and spatially heterogeneous scenes such as tropical forests, 3D models are more suitable since they account for vertical and horizontal heterogeneity, and include shading and multiple scattering [46]. 3D RTMs reconstruct the architecture of a medium through geometric primitives such as triangles, cubes, cones, spheres or cylinders [47]. The combination of these elements with the other inputs of the RTMs allows generating 3D scenes of different levels of complexity.
Simple geometric structures derived from dendrometric measurements were initially used to approximate the forest structure in the first studies using 3D RTM [47]. The recent development of light detection and ranging (LiDAR) provided an unprecedented level of detail to characterize forest 3D architecture [48]. Schneider et al. [49] evidenced the interest of using detailed 3D structure of vegetation derived from LiDAR compared to simpler 3D representations of individual trees defined as geometric shapes and located in the scene based on field inventories: they obtained increased accuracy when comparing experimental imaging spectroscopy with simulations using LiDAR-derived 3D structure.
In the current era of deep learning, RTMs are also strong candidates in the perspectives of data augmentation [50] to alleviate challenges related to data sparsity in order to understand interactions between vegetation and solar radiation, to test potential and limits of methods, and to prepare future missions. DART (Discrete Anisotropic Radiative Transfer) [51,52] is currently one of the most comprehensive 3D RTM for simulating remotely sensed data, in particular in the context of forest studies [42,53]. However, the proper parameterization of DART in order to account for forest structural and radiometric properties requires particular care, and needs to be tuned for the application considered. Malenovský et al. [54] showed the influence of including woody elements in 3D simulations in order to retrieve Leaf Area Index (LAI) and simulate canopy reflectance for a Norway spruce forest. Strong assumptions on the homogeneity of LOP within species or individuals are also usually made when using 3D models to study heterogeneous vegetation. Schneider et al. [49] simulated imaging spectroscopy with DART on a temperate mixed forest in order to study the potential of upscaling leaf level observations, and put efforts on characterization of the vertical heterogeneity of LOP for each tree species, but used the same LOP for a given layer and all individual trees from the same species. Ferreira et al. [55] studied the contribution of DART simulations to the retrieving of structural and chemical properties of individual tree crowns (ITCs) in dry tropical forests, and assumed homogeneity of LOP in ITCs. The assumptions made by the authors were appropriate to meet their objective. However, accounting for intra-species and intra-crown variability of LOP is potentially crucial for the production of realistic simulations aiming at exploring the potential of methods dedicated to biodiversity monitoring since there is spectral variability within crowns, within species and among species in tropical forest [56]. Thus, the hypothesis of uniform LOP at either crown or species scale may lead to dramatic underestimation of spectral heterogeneity in simulation of forest stand scale remotely sensed data, making these simulations unsuitable for applications built upon the fine analysis of spatial heterogeneity of spectral information, such as biodiversity mapping.
The lack of information on the non-photosynthetic vegetation fraction NPVf and an inadequate characterization of the diversity of LOP have been shown to contribute to the relatively low radiometric agreement between simulation and experimental data [49,54,57]. Simulating the variability in canopy reflectance is challenging in many ways, as it requires definition of fine scale 3D structure and heterogeneity of LOP. While Airborne Laser Scanning (ALS) potentially fulfills the requirements for the definition of the 3D structure, the definition of LOP also needs to be addressed using indirect estimation, as sampling within and among species variability in optical properties in hyperdiverse ecosystems is not feasible. To address those challenges, we designed a methodology aiming at simulating spatial heterogeneity in spectral data over a complex tropical forest using the DART model. We used ALS data in order to reconstruct 3D forest architecture.
Then, we defined scenarios to account for NPVf using different proxies based on 3D structure or canopy reflectance. Once an optimal scenario to characterize NPVf was defined, we performed simulations accounting for LOP variability at different scales (pixel, crown, species). Finally, we compared simulations to airborne imaging spectroscopy using statistical metrics, dissimilarity metric and supervised classification methods. The main objective of this study was to define a framework for producing realistic simulations in the perspective of supporting methodological development for biodiversity assessment using 3D modeling, considering: i) the effects of incorporation of NPVf into simulation and ii) the influence of LOP variability on the spectral response simulated by the DART model.

Study area
The study area is located in the coastal part of French Guiana at Paracou research station (5°18′N, 52°53′W). It is characterized by a lowland terra firme rainforest. The forest is dominated by Leguminosae-Caesalpinioideae, Lecythidaceae, Chrysobalanaceae, and Sapotaceae families. Field plot network inventories carried out almost every year for more than 35 years highlight a large specific diversity with more than 750 tree species identified over 118.75 ha. A complete description of the Paracou experimental site can be found on https://paracou.cirad.fr/website.

Remotely sensed data
In September 2016, an airborne campaign funded by the Centre National d'Etudes Spatiales (CNES) and the Center for the study of biodiversity in Amazonia (labex CEBA) led to the acquisition of imaging spectroscopy, LiDAR data, and very high spatial resolution RGB imagery. The acquisitions were performed in clear sky conditions at an average flying altitude of 890 m.

Imaging spectroscopy
Imaging spectroscopy data was acquired with a Hyspex VNIR-1600 sensor (Hyspex NEO, Skedsmokorset, Norway) over 160 spectral bands in the visible and near infrared (VNIR) regions, ranging from 415 nm to 993 nm. Spectral information corresponding to spectral bands > 902 nm were discarded because of the low signal-to-noise ratio. Images were orthorectified and geo-referenced using the PARGE software (http://www.rese.ch/products/parge/) with the canopy DSM (Digital Surface Model) produced from the LiDAR point cloud. Atmospheric corrections were performed using ATCOR-4 software [58] to produce Top of Canopy (TOC) reflectance with 1-m spatial resolution. Airborne acquisition over the full site and its surroundings resulted in 23 flight lines covering 14 Km 2 .

Airborne laser scanning (ALS)
ALS data was acquired with a RIEGL LMS Q780 simultaneously with imaging spectroscopy data. The sensor operated with a frequency of 400 kHz and a scanning angle of ± 30° which ensured an average pulse density of 33 pulses/m². A point cloud filtering was applied to remove non-valid points.

RGB imagery
Very high spatial resolution (VHSR) imagery was acquired during the same airborne campaign using an iXu 180 Phase One camera. This imagery includes red, green and blue channels and the final spatial resolution of the orthorectified mosaic was 10 cm.

Ground information
2.3.1. Field plot network and selected sites for simulation The Paracou site hosts an experimental setup including multiple permanent plots contributing to long-term monitoring of the influence of human induced disturbances since 1986 [59]. Species inventories combined with accurate tree trunk location have been performed every 1-2 years for each plot for more than 35 years for oldest plots.
We selected two subplots in the control undisturbed forest plots, dimensioned to fulfill the different aspects of our analysis: • A plot of 80 x 80 m, hereafter referred to as Site A was dedicated to the adjustment of statistical models for the estimation of vegetation biochemical properties; • A plot of 300 x 300 m, hereafter referred to as Site B, was dedicated to the comparison of the simulations with experimental data.
The extent of the airborne imaging spectroscopy acquisition and the location of the plot network as well as Sites A and B are presented in Figure 1.

Individual tree crown delineation
Manual delineation of emerging crowns identified from the inventories was conducted, based on visual interpretation of the VHSR image, combined with a canopy height model derived from ALS. This delineation was validated during field surveys [60][61][62]. More than 2500 ITCs were delineated in the area covered by the imaging spectroscopy acquisition, but our analysis focused on ITCs used as support to the modeling activities performed on site B. We selected the nine most abundant species from these delineated crowns in site B, corresponding to 162 ITCs and 10 262 pixels. Table 1 provides basic statistics on size and distribution of ITCs per species.

Overview of the methodology
Our main objective is to define a simulation framework aiming at simulating realistic remote sensing data measured with optical sensors over tropical forests. To reach this goal, we focused on identifying the optimal method for the introduction of two major contributors of the heterogeneity measured with remote sensing, NPVf and LOP. Figure 2 provides the general workflow of our study. During the first step, the global parameterization of DART simulations was defined for all the factors remaining unchanged for all simulations, based on experimental data. This global parameterization included sensor characteristics, geometry of acquisition, topography definition as well as 3D forest reconstruction from ALS data. The second step aimed at building different scenarios for the integration of two factors into the simulations, NPVf and spatial variability in leaf optical properties. NPVf was based either on the introduction of woody elements (Wood Area Density, WAD) into the 3D reconstruction, or on the addition of brown pigments (Cbrown) into leaf constituents. The spatial variability in leaf optical properties was defined at different levels of organization (pixel, crown, species) based on inventoried and delineated tree crowns. The third step consisted in running DART model in order to produce simulated imaging spectroscopy data for each scenario. Finally, the fourth step aimed at comparing simulations with airborne acquisitions at different levels: we compared experimental airborne acquisitions with simulations based on pixel reflectance, spectral dissimilarity within and among species, and species discrimination performed with a supervised approach.

Leaf optical properties modeling with PROSPECT
We used the leaf model PROSPECT [63,64] to simulate leaf optical properties. PRO-SPECT is a leaf model simulating directional-hemispherical reflectance using biophysical and biochemical parameters. It is based on Allen's generalized plate model [65]: a leaf is modeled as a pile of identical elementary layers characterized by an absorption coefficient and a refractive index defined over the optical domain from 400 nm to 2500 nm. Here, we used PROSPECT-D [63] including three main pigments that control the absorptive properties of fresh leaves in the visible (VIS) domain, i.e., chlorophylls a+b (Cab), carotenoids (xanthophylls and carotenes) (Cxc), and anthocyanins (Canth). Additional constituents accounted for by PROSPECT-D include brown pigments (Cbrown), which is used to simulate senescent leaves, and contributes to leaf absorption in the VNIR, leaf mass per area (LMA) and equivalent water thickness (EWT), which are mainly contributing to leaf absorption in the short wave infrared domain. The scattering properties of the leaf are controlled by the leaf structural parameter (N), which corresponds to the number of elementary layers in the leaf, separated by N-1 air layers.

3D modeling of Canopy reflectance with DART
The DART (Discrete Anisotropic Radiative Transfer) model [51,52] was used to simulate airborne imaging spectroscopy. DART simulates the radiative transfer in the 3D "Earth-Atmosphere" system from the visible to thermal infrared domain of the electromagnetic spectrum. It is used to study both natural and urban landscapes, and allows for the integration of atmospheric and topographic effects. The DART scenes were based on 3D representations made of individual volume elements (voxels), each of them characterized by a set of properties such as vegetation density, leaf angle distribution, optical properties corresponding to leaves, trunk and branches, soil and any relevant elements included in the scene. Simulations were performed with the package pytools4dart (https://pytools4dart.gitlab.io/pytools4dart/) a python API aiming at preparing, running and analyzing DART simulations, specially dedicated to the design of large and complex simulations.
DART was set in flux tracking mode to simulate TOC reflectance, with the sun and the atmosphere as the only radiation sources. Topography was defined with a digital surface model produced from ALS point clouds. A generic litter spectrum was used for the optical properties of the soil/background. Simulations were performed for 135 bands ranging from 415 nm to 902 nm, following instrumental characteristics of the airborne acquisition for sites A and B. Sensor height and field of view, solar zenith and azimuth angles as well as spatial resolution were set according to the geometry of acquisition during the airborne campaign.

Defining scenarios for the integration of LOP heterogeneity in 3D modeling
Spectral variability in canopy reflectance among species, within species, and within individual trees has been documented [56] and multiple factors such as position in canopy or topography contribute to this variability. In our study, we specifically focused on the influence of spatial variations in LOP on this spectral variability: we defined three scenarios for the integration of LOP in radiative transfer simulations at pixel level (LOPpixel), crown level (LOPITC), or species level (LOPspecies). Then we compared simulated canopy reflectance with corresponding experimental acquisitions in order to describe the influence of LOP variability on spectral heterogeneity of canopy reflectance and assess the degree of realism of the simulation.
The definition of LOP at different scales could not be performed using field sampling for logistic reasons. Instead, we implemented a method to estimate a set of relevant chemical constituents from canopy reflectance: here we focused on leaf pigments, and estimated the combination of pigment content for each pixel. Then we used the pigment content estimation at pixel level, or its average at crown or species level using the species inventory and the crown delineation, to simulate and assign LOP values in the simulation scene depending on the scenario.
Following the parameterization of 3D structure and leaf chemistry defined, simulations were performed in order to produce imaging spectroscopy data for the different scenarios. Finally, a set of indicators was computed in order to compare the simulated scenarios with the experimental imaging spectroscopy data, and assess the realism for each of them.

Definition of 3D structure of the canopy
The 3D forest structure was derived from ALS acquisition: the ALS point cloud and associated trajectory file were voxelized using the AMAPvox v.1.7.3-BETA software (http://amap-dev.cirad.fr/projects/amapvox). AMAPvox tracks every laser pulse through a 3D grid (voxelized space) from the top of the canopy to the last recorded hit. The effective sampling area of each laser pulse (or fraction of pulse in case of multiple hits) is computed from the theoretical beam section (a function of distance from laser and divergence of laser beam) and the remaining beam fraction entering a voxel. In case more than one hit is recorded for a given pulse then it is assumed that the surface intercepted by each hit (the "hit area") is equal to the inverse of the total number of hits times the beam section. This information is combined with the optical path length of each pulse entering a voxel to compute the local transmittance per voxel. Different estimation procedures are provided in the AMAPVox software. In the present study we used a numerical resolution to compute for each voxel the transmittance -or gap probability Pgap, from the equation below: With the optical potential path length of pulse i (from entry point to the expected exit point) the cross section of pulse i, computed at voxel center and being the beam surface fraction entering and outgoing the voxel The equation above expresses that the estimated per unit length of a voxel is the value that yields the light transmission rate (given the actual path length and effective surface of all contributing pulses) equal to the observed transmission rate. The gap fraction probability per unit path length per view direction depends only on the foliage projection ratio and the plant area density PAD (m 2 .m -3 ).
G(θ)projection function will depend on the leaf inclination distribution function (LIDF) which is also an important factor among canopy structure properties influencing canopy reflectance [48,49,66]. However available data did not allow to include estimated LIDF in our simulations, and we used a unique LIDF corresponding to spherical distribution, the most widely used since it offers straightforward computation of the projection function, equal to 0.5 in any direction [66].
Here, we used a uniform 1 x 1 x 1 m voxel size matching with the spatial resolution and pixel grid of the imaging spectroscopy acquisition. Some post-processing steps are required to adjust the raw 3D distribution of PAD obtained from AMAPvox [67][68][69][70] and ensure accurate description of horizontal and vertical structure. These adjustments were performed as follows: • A hierarchical re-estimation of the transmittance based on a mixed linear model was performed in order to minimize the uncertainty associated with poor sampling of lower canopy voxels resulting from the gradual extinction of the laser beam (occlusion) [70].
• A reduction factor of 0.8 was applied to the re-estimated PAD in order to compensate for the limited penetration achieved by the LMSQ780 lidar when flying at 900m. This factor was set to fit the extinction profile obtained over the same area while operating at lower altitude (450m) (G. Vincent, in prep).

Integration of NPVf into simulations
Each voxel in the 3D scene is characterized by a specific PAD. PAD is an important factor driving the interactions of electromagnetic radiation with vegetation [36,71], which combines the contribution of various elements including leaves of all development stages (from juvenile to senescent and dead leaves), branches and stems. The vertical integration of the PAD corresponds to the Plant Area Index (PAI). PAI is defined as the sum of Leaf Area Index (LAI) and Wood Area Index (WAI), under the assumption of homogeneous canopy, excluding clumping effects [72]. LAI is a key parameter to consider for physiological, ecological and climatological studies [73]. Old-growth forests are characterized by a non-negligible fraction of NPVf corresponding to dead leaves and woody components like standing deadwood, twigs and branches distributed within tree crowns [74]. Zhu et al. [75] showed an overestimation of LAI estimated by terrestrial LiDAR of 3.0 to 31.9 %, caused by woody elements, with differences observed among forest types. Calders et al. [76] analyzed the relationship between effective PAI (ePAI), effective LAI (eLAI), and effective WAI (eWAI), which corrects for clumping effect, using the librat Monte Carlo ray tracing model. They obtained a linear relationship between eLAI and ePAI-eWAI. However, no strict complementarity was obtained between these three area indices, as woody structure is preferentially obscured by the leaves. Accounting for the relative effect of these woody elements is particularly challenging when using radiative transfer modeling [77] since they are important photon absorbing and scattering components [36], hence both their spectral response and their 3D distribution in a forest scene should be considered in radiative transfer approaches.
In our study, we compared two approaches to take into account NPVf in our simulations, based on different hypotheses: the first approach is based on the disaggregation of PAD into LAD and WAD, while the second approach assumes equality between PAD and LAD and uses brown pigment content in leaves to account for the influence of NPVf on canopy reflectance. The former will be referred to as NPVfWAD, while the latter will be referred to as NPVfCbrown. We also compared these two scenarios accounting for NPVf with the case where PAD is exclusively defined as leaf materials without brown pigments, hereafter referred to as NPVf0.

Estimation of NPVfWAD
The definition of WAI from the voxelized mock-up is based on two hypotheses. First, PAD is assumed to be the sum of Wood Area Density (WAD) and LAD: although literature showed that this hypothesis is incorrect when accounting for clumping [72,78], we used 3D forest representation based on individual voxels characterized by turbid medium with no clumping. Second, the contribution of WAD to PAD is assumed to be proportionally constant. This assumption is not validated by experimental observations, but this simplification is imposed by the impossibility to accurately differentiate wood and leaves from ALS data.
In order to adjust the proportion of WAD in each voxel, we first identified a set of ITC corresponding to defoliated and foliated individuals based on visual analysis of VHSR imagery. This visual identification resulted in 10 ITCs unambiguously identified as defoliated, and 20 ITCs unambiguously identified as fully foliated and non-senescent. The visual identification was not assisted by inventory data, so our selection included a random selection of species. Then, we extracted the PAD corresponding to the three uppermost voxels from the 3D mockup, and we computed the mean PAD from foliated and defoliated voxels. We defined the reference NPVfWAD as the ratio between the mean PAD from defoliated voxels and the mean PAD from foliated voxels, assuming that in average the WAD of the defoliated crowns is similar to the WAD of foliated crowns. In order to account for woody structure being preferentially obscured by the leaves, we also tested an intermediate scenario in which NPVf is half of NPVfWAD. This intermediate scenario will be referred to as NPVfWAD/2.

Estimation of leaf brown pigment content
In the second approach, we assumed equality between PAD and LAD, and used a spectral index in order to estimate the brown pigment content (Cbrown) corresponding to each pixel. This approach is based on the hypothesis that the radiometric contribution of non-photosynthetic elements in a voxel can be compared to the radiometric contribution of Cbrown, which absorbs in the VNIR domain. Brown pigments are characterized by a broad absorption domain, with maximum absorption in lower wavelengths of the VIS domain, gently decreasing until 1100 nm [79]. The introduction of brown pigments results in smoothing of the green reflectance peak, and more remarkably a less abrupt increase in reflectance in the red-edge. The influence of Cbrown in the red-edge domain can be used to define a spectral index linking canopy reflectance to leaf brown pigment content. The definition of this spectral index and the procedure used to adjust the regression model linking canopy reflectance to leaf pigment content is described in the next section.

Simulation of leaf optical properties through the estimation of leaf pigment content
Leaf pigment content controls LOP and canopy reflectance in the VIS region of the electromagnetic spectrum as well as in the red-edge domain and first part of the NIR domain [36]. Here, our first goal was to define a strategy to assign pigment content to each pixel, crown or species in our simulations, based on canopy reflectance. A variety of methods can be used to assess foliar pigment content from remotely sensed canopy reflectance acquired with imaging spectroscopy. Verrelst et al. [80] categorized these methods into four categories: (i) parametric regression, including vegetation indices, shape indices and spectral transformations, and assuming a simple relationship between spectral data and vegetation biophysical variable, such as a linear monovariate relationship; (ii) nonparametric regression, based on data-driven techniques that directly link spectral data and biophysical variables (e.g. linear and nonlinear machine learning regression algorithms); (iii) inversion of RTM using numerical optimization and look-up table approaches; and (iv) hybrid methods, which combine RTM for the simulation of a training dataset with parametric or nonparametric regression methods.
Parametric regressions based on vegetation indices and linear or polynomial models are the most straightforward means of estimating pigment content [81][82][83]. However, the adjustment of such empirical models requires a large amount of field data. Hybrid approaches provide a means to overcome this dependency on field data: they are based on the simulation of a training dataset (corresponding to LOP or canopy reflectance) using appropriate physical models and sampling strategy, in order to adjust a regression model applicable to experimental data. Therefore, we chose a hybrid approach in this study, in order to train regression models linking spectral indices with leaf pigment contents, including Cab, Cxc, Canth and Cbrown. This approach is based on two steps: First, the regression models were adjusted based on a simulation including a broad range of pigment content. Here, we adjusted a linear model specific to each pigment and each NPVf scenario, based on simulations performed on Site A. Then the linear models were applied on experimental data corresponding to Site B, in order to define specific pigment content for each NPVf scenario.
3.6.1. Adjustment of models for the estimation of leaf pigments (Site A) We used three spectral indices to estimate Cab, Cxc and Canth from remotely sensed data: The Chlorophyll Index (CI) (Equation 3), the Carotenoids Reflectance Index (CRI) (Equation 4) and the Anthocyanin Index (ARI) (Equation 5) [84]. These spectral indices were applied on imaging spectroscopy for various ecosystems including forest [28,29,85]. Due to the lack of documented spectral index for Cbrown in the literature, we defined a Brown Pigment Reflectance Index (BPI) (Equation 6) based on reflectance measured in the rededge and NIR domains. The spectral band selected in the red-edge was centered on 750 nm, which is characterized by negligible absorption of Cab compared to Cbrown.
Where Rλ1-λ2 corresponds to the mean reflectance in the spectral domain defined by λ1 and λ2.
DART was used to simulate the hyperspectral acquisitions on site A coupled with PROSPECT-D to simulate LOP. We followed the global parameterization setting for instrumental configuration, the structure of the 3D mockup was parameterized as described in Section 3.4. A sampling strategy based on experimental distribution and co-distribution of leaf chemistry [86] was designed in order to define a combination of leaf pigments to be associated with each pixel, i.e. vertical column of voxels, of the 3D mockup. We then generated 6400 combinations of Cab, Cxc, Canth and Cbrown following a normal multivariate distribution for Cab and Cxc, a log normal distribution for Canth, and a normal distribution for Cbrown. Each of these 6400 combinations was assigned to a pixel from the simulation corresponding to site A.
The other input parameters of PROSPECT-D, N, EWT and LMA, were set to a constant value. The leaf structure parameter N influences leaf scattering properties and the ratio between reflectance and transmittance, but does not influence leaf absorption. Hence, the N parameter was set to 1.8, assuming that the vegetation scattering induced by the detailed canopy structure would provide sufficient variability. The influence of EWT and LMA on canopy reflectance is limited in the VNIR spectral domain, thus they were set to 12.9 mg.cm -2 and 12 mg.cm -2 respectively, which correspond to the average values reported by Féret et al. [86] for a leaf dataset collected over the Paracou research station.
The linear regression models between leaf pigment content and the spectral indices (computed from simulated canopy reflectances) were trained using sunlit pixels only, i.e. pixels with more than 20 % reflectance in the NIR spectral band centered on 800 nm, in order to reduce confusions induced by low illumination in further comparisons [23].

Application of models for the estimation of leaf pigments (Site B)
The regression models adjusted on site A were applied to the pixel of the airborne hyperspectral observations of Site B in order to estimate pigment content site B 3D mockup. Same as for site A, all voxels in the same vertical columns shared the same pigment content. The pixel-scale pigment content was then aggregated at crown scale following delineated ITCs and at species scale using the taxonomic identification of these ITCs, using the average value of sunlit pixels. A generic spectrum simulated using a combination of canopy pigment was used to characterize pixels of non-delineated ITCs in Species and ITC scenarios. Finally, for each level of LOP heterogeneity (pixel, crown, species), the four NPVf scenarios (NPVf0, NPVfWAD, NPVfWAD/2, NPVfCbrown) were produced, resulting in 12 simulations of imaging spectroscopy data corresponding to Site B.

Comparing simulations with airborne acquisitions
We compared airborne acquisitions with simulated imaging spectroscopy data based on different levels of information. It includes radiometric comparison between measured and estimated canopy reflectance, as well as comparison of crown/species scale spectral variability between measured and estimated canopy reflectance, species discrimination using supervised classification algorithms and spectral dissimilarity within ITCs, between ITCs of the same species and between ITCs of different species using spectral angle metric.

Radiometric comparison
Different statistical criteria were used to compare the observed ( , ) and simulated TOC reflectance ( , ) in the VNIR spectral domain. These include the root mean square error (RMSE) (Equation 7) and the coefficient of determination (R 2 ) (Equation 8), computed as follows for each spectral band of the sensor:

Analysis of spectral dissimilarity
Tropical tree species discrimination using supervised methods remains challenging, and assuming a direct link with within-and among-species variability of canopy reflectance is a simplification [87]. A number of statistical metrics of variability such as the coefficient of variation, the spectral angle or the minimum distance from the spectral centroid can be used to relate spectral diversity to plant diversity [22]. In this study, the spectral angle [88], defined in (Equation 7), was used to compare spectral dissimilarity between experimental data and simulated data. The spectral angle was computed on sunlit pixels at three level using the 135 spectral bands: (i) within crown for all pairwise combinations of pixels of a given crown, (ii) within species for all pairwise combinations of crowns of a given species, and (iii) among species, for all pairwise combinations of species. The comparison between the dissimilarity matrices of observations and simulations was done at each level with Pearson's correlation coefficient and RMSE.
Where Si and Sj are the spectra of element i and j, and λa and λb are the wavelength boundaries i.e., 415 to 902 nm in this study.

Species discrimination
A large number of publications evidenced the interest of imaging spectroscopy data for the discrimination among species, even for tropical forests showing strong diversity [24,61,87,[89][90][91][92]]. On a higher ecological level, this capacity to discriminate among species based on spectral information is also a fundamental property for methods aiming at estimating diversity indices based on the heterogeneity of spectral information [28,29,34]. In this study, we used the performances of species discrimination based on supervised classification as an indicator of the realism of our simulated imaging spectroscopy data. Species discrimination applied in real situations inherently leads to successful classification for species characterized by high discriminability, and confusion between species showing lower discriminability. This capacity to accurately identify a species or misidentify it as another species is explained by the expression of vegetation traits. In our study, we applied supervised classification for species identification on both experimental and simulated imaging spectroscopy data, hypothesizing that realistic simulations would lead to similarity in the confusion matrix derived from classification when compared to the confusion matrix derived from experimental data. We used the linear discriminant analysis (LDA) as a supervised classification algorithm in order to perform species classification [61,87,90].
LDA is a parametric classification algorithm based on the hypothesis of multivariate normality of input variables, and aims to maximize the between-class variance and minimize the within-class variance, through a linear discriminant function. In our study, we performed classification on ITCs corresponding to the nine most abundant species detailed in Table 1. For each species the ITCs were randomly distributed between two folds of equal size: one for training and the other for test. The split is made at the crown level instead of the pixel level to avoid that pixels of the same crown, thus with a possible spatial correlation, could be found in both training and test folds, which would wrongly increase the test scores. This resulted in 4 to 23 ITCs/species used to train the LDA model, and 4 to 22 ITCs/species used to test the model. The procedure was applied on imaging spectroscopy data, with training and test datasets built from the same pixel selection for experimental and simulated data, and repeated 100 times in order to account for the influence of the random selection. The Pearson's correlation coefficient was used to compare the average confusion matrices of tests on observations against simulations. Two correlation analyses were performed, one to compare the success of classification using the diagonal elements of the confusion matrix, and another one comparing confusion rate between species and using off-diagonal elements.

Definition of the 3D structure of the forest
The PAI was computed for each pixel in site B, first from the initial 3D mock-up produced with AMAPVox, after application of the hierarchical re-estimation of transmittance, and finally after application of the attenuation factor. The mean PAI following was 9.71 m²/m² for Site A and 9.4 m²/m² for Site B.

NPVf derived from WAD
The PAD of three uppermost voxels of the selection of foliated and defoliated ITCs was averaged in order to analyze their distribution. The mean PAD corresponding to defoliated and foliated tree crowns was 0.33 m 2 .m -3 and 1.54 m 2 .m -3 , respectively ( Figure 3). NPVfWAD, defined as the ratio of these two mean densities, was then set to 21% in simulations. Other studies reported NPVf ranging between 10% and 50% depending on forest type and stage [55,80,93]. In order to account for potential shadowing of woody elements by leaves [61], we also tested a scenario corresponding to 50% of NPVfWAD hereafter referred to as NPVfWAD/2. For each of these scenarios, the optical properties corresponding to generic wood optical properties were assigned to a constant proportion of PAD defined as NPVfWAD or NPVfWAD/2. The complementary fraction of PAD corresponding to LAD was characterized by LOP with pigment content defined as described in Section 4.3.

NPVf derived from leaf brown pigment content
A linear regression model was adjusted in order to describe the relationship between BPI and Cbrown based on DART simulation over Site A. The Pearson's correlation coefficient (r) corresponding to this linear relationship was 0.67 (Figure 4a). The linear model was then applied to sunlit pixels from Site B in order to account for non-photosynthetic vegetation in the scenario considering the PAD as equivalent to LAD. The value of Cbrown ranged between 0 and 0.6 arbitrary units (Figure 4b). The scenario accounting for NPVf through Cbrown is referred to as NPVfCbrown hereafter.

Impact on the adjustment of the relationship between spectral indices and pigment content
Once DART simulations were performed on Site A for the different scenarios accounting for NPVf, the relationships between pigments content and their corresponding spectral indices were studied ( Figure 5). Overall, a linear relationship with a positive correlation was observed between pigment content and corresponding spectral index, for all pigments and all scenarios accounting for NPVf. The relationship between Cab and CI showed the strongest correlation, with Pearson correlation coefficient ranging between 0.89 and 0.93. However, the correspondence between Cab and CI strongly varied depending on the scenario used to account for NPVf: NPVfWAD and NPVfWAD/2 resulted in dramatic decrease in the range of CI ( Figure 5, X-axis), which is explained by the decreased effect of chlorophyll absorption on the signal when mixing photosynthetic elements with nonphotosynthetic elements.
The relationship between Cxc and CRI showed lower correlation, with Pearson's correlation coefficient ranging between 0.55 and 0.65, with maximum correlation obtained for NPVf0, and comparable correlation for the three other scenarios. Like for CI, the range of variation observed for CRI decreased with NPVfWAD and NPVfWAD/2 due to reduced contribution of photosynthetic vegetation.
Finally, the relationship between Canth and ARI showed relatively different behavior: The Pearson's correlation coefficient was similar for NPVf0 and NPVfCbrown (Pearson's r= 0.83 -0.82), while it increased to 0.91-0.90 when including wood optical properties in the voxels, along with a reduced range of variation of ARI, as for the two other indices.
While all spectral indices appeared to be sensitive to the pigment content they were initially designed to be used with, the linear models adjusted based on DART simulations are expected to result in very different pigment content when applied to experimental data, depending on the method used to account for NPVf.

Implication of the application of the models on experimental data
The regression models adjusted based on DART simulations performed on Site A ( Figure 6) were applied on sunlit pixels from Site B in order to estimate pigment content for each NPVf scenario. The distribution obtained for Cab and Cxc was very similar for NPVf0 and NPVfCbrown. Cab ranged from 2.9 to 82.7 µg.cm -2 with mean of 35.2 ± 10.1 µg.cm -2 for NPVf0 and from 0.9 to 93.7 µg.cm -2 with a mean value of 38.5 ± 11.7 µg.cm -2 for NPVfCbrown. The two scenarios NPVfWAD and NPVfWAD/2 resulted in increase in Cab ranging from 1.6 to 105.4 µg.cm -2 with a mean of 43.7 ± 13.1 µg.cm -2 , and from 1.1 to 133.6 µg.cm -2 with a mean value of 54.8 ± 16.7 µg.cm -2 , respectively.
The influence of the different NPVf scenarios on the estimation of Cxc was comparable to the results obtained with Cab. Cxc ranged from 1.5 to 17.9 µg.cm -2 with a mean value of 7 ± 2.4 µg.cm -2 for NPVf0 and from 1.7 to 23.7 µg.cm -2 with a mean value of 9.1 ± 3.2 µg.cm -2 for NPVfCbrown. The distribution of Cxc strongly increased for NPVfWAD and NPVfWAD/2 scenarios.
The influence of the different NPVf scenarios on the estimation of Canth differed from the influence observed for photosynthetic pigments. Canth distributions including negative values were obtained for all scenarios with various occurrences, suggesting discrepancies between simulated and measured reflectance for the spectral domain used to compute ARI. These negative values were set to zero when introduced into DART simulations.

Influence of NPVf on simulated reflectance
Imaging spectroscopy and DART simulations were compared using spectral R 2 and the RMSE for each NPVf scenario (Figure 7). We assumed that higher R 2 and lower RMSE would indicate higher radiometric realism. Overall, NPVfCbrown and NPVf0 scenarios resulted in similar R 2 over all the spectral domain, with stronger correlations than those obtained for NPVfWAD and NPVfWAD/2 over the full spectral domain. The scenario corresponding to NPVfWAD resulted in lowest R 2 overall. The spectral profile of R 2 was similar among NPVf scenarios: lower correlations were obtained in the far-blue (<450 nm), green (520 nm), red-edge and NIR domains, while higher correlations were obtained between blue and green domains (500 nm), and in the yellow-red domains (570 nm -680 nm).
The choice of NPVf scenario resulted in strong differences in RMSE between measured and simulated reflectance in the red-edge/NIR domains. The NPVfCbrown scenario resulted in lower RMSE compared to other scenarios in the full spectral domain, except in the red domain. Spectral bands located in the blue, green, red, red edge and NIR spectral domains were selected in order to analyze the influence of NPVf on the simulated reflectance ( Figure 8). The influence of NPVf on reflectance was particularly strong at 491 nm (blue) and 666 nm (red), corresponding to domains of absorption of photosynthetic vegetation. In these regions, an underestimation of reflectance was observed for NPVf0. The integration of NPVf using WAD led to an increase of simulated reflectance up to an overestimation for NPVfWAD. For NPVfCbrown, reflectance was underestimated at 491 nm and 666 nm. In the green domain (560 nm), the effect of NPVf was moderate, except for NPVfCbrown which removes the positive bias observed with the other scenarios. An overestimation of reflectance was observed for NPVf0 at 720 nm (red edge) compared to scenarios where NPVf were considered. Simulated reflectance was overestimated at 800 nm overall but no significant effect of NPVf was detected. Given these results, we decided to only use the NPVfCbrown in the following analyses.

Influence of LOP variability on simulated reflectance
Spectral R 2 and RMSE were used to compare reflectance from airborne acquisitions and DART simulations for each spectral band and each LOP level, when selecting NPVfCbrown. LOPpixel outperformed LOPITC and LOPspecies in the visible region based on R 2 and in the NIR based on the RMSE (Figure 9). Relatively similar results were obtained for other NPVf scenarios (results not shown).

Analysis of spectral dissimilarity
The spectral dissimilarity was computed within crowns, among crowns corresponding to the same species and among crowns corresponding to different species using the spectral angle. DART simulations were compared to imaging spectroscopy to test if the spectral dissimilarity observed in imaging spectroscopy data was correctly simulated for the different LOP scenarios ( Figure 10 and Table 2). When defining LOP per species in DART (LOPspecies), low agreement was obtained between spectral dissimilarity computed from imaging spectroscopy and from simulations for all levels (within crowns, within species and among species), with systematic underestimation. Intermediate performances were obtained when defining LOP per ITC in DART simulations. LOPITC scenario showed a systematic underestimation of spectral dissimilarity within crowns and a good agreement between spectral dissimilarities within species and among species. The definition of LOP per pixel in DART simulations outperformed other scenarios, with strong agreement of spectral dissimilarity when comparing experimental data with simulations for all levels.

Assessing the potential of species discrimination
Performances for species classification based on LDA were compared when applied on airborne imaging spectroscopy and DART simulations ( Figure 11 and Table 3). The simulation corresponding to LOPspecies resulted in an overestimation of the potential of species discrimination compared to airborne imaging spectroscopy, with close to 100% overall classification accuracy. The confusion matrix derived from classification applied on LOPITC and LOPpixel scenarios showed stronger consistency with the confusion matrix derived from classification applied on experimental. Overall, LOPpixel outperformed LOPITC, with particular ability to reproduce confusion among species.  Table 3. Pearson's correlation coefficient comparing LDA classification of DART simulations to imaging spectroscopy for LOP level scenarios. The consistency in correct classification and misclassification is compared between classification results obtained from airborne acquisition and classification results obtained from simulations.

LOPspecies LOPITC LOPpixel
Correct classification 0.05 0.86 0.92 Misclassification -0.04 0.71 0.79 Figure 12 allows for visual comparison between imaging spectroscopy and the DART simulation of site B with the LOPpixel/NPVfCbrown scenario. Similar patterns and object dimensions were observed after stretching the nadir true color composites of DART simulation ( Figure 12b) and imaging spectroscopy (Figure 12e). The nadir true color composites ( Figure 12a) and imaging spectroscopy (Figure 12d) highlight the need for spectral characterization of defoliated crowns. This aspect was not considered in this study. In addition, the comparison in Figure 12f highlights the need for accurate characterization of biophysical and structural parameters which particularly affect the NIR region. This, in addition to the improvement of pigment contents estimates, will potentially increase radiometric agreement between simulations and imaging spectroscopy.

3D forest reconstruction
The 3D forest structure was reconstructed by discretizing ALS data in 3D space using AMAPVox software. By applying Beer-Lambert's turbid medium approximation, local transmittance enables to estimate PAD and PAI. Infravoxel clumping may have caused an underestimation of LAI and PAI [75,94]. Qu et al. [95] estimated the PAI of a moist tropical forest at Fazenda Cauaxi in the Paragominas Municipality of Pará State, Brazil. They reported a mean PAI of 6.94 m²/m². Another study was conducted by Pimmasarn et al. [67] in Khao Yai National Park in central Thailand. They obtained a mean PAI of 8.02 m²/m². These results are slightly lower than the PAI reported here (mean PAI of 9.4 m²/m²). On the opposite, Vincent et al. [70] mapped PAI at the same site (Paracou experimental station) by ALS (using a different laser system, i.e. a Riegl LMSQ560) and reported a mean PAI of 13.2 m²/m². In that study the authors also compared simulated and measured light levels below canopy which revealed a slight underestimation of the ALS estimated transmittance [70]. Moreover, it is worth noting that the formula used to compute voxel transmittance from ray tracing in Vincent et al. [70] (equation 7) can be shown to be biased by the fact that it considered an average optical path length rather than the actual path lengths as was done here. The bias might be responsible for the underestimation in transmittance (and overestimation of PAI) observed at that time. This may indicate that the present estimate is more accurate than the one previously published for the site.
An independent measurement of PAI on site A (80 x 80 m) was conducted using high density TLS in 2016 (G. Vincent and F. Heuschmidt (2018) unpublished). The plot was voxelized at 50cm resolution using AMAPvox. The estimated plot level PAI was found to be equal to 11.4 m²/m².
Finally, it should be noted that the uncertainty in PAI which to a large extent is related to the limited sampling density in the lower canopy due to occlusion. The upper canopy is much better sampled and the PAD of the top canopy layers which indeed are of prime interest here are probably more accurately estimated than PAI (PAD integrated over the entire profile).

Influence of NPVf on leaf chemistry estimate
A set of linear models was adjusted to estimate pigment contents from spectral indices, for each of the scenarios accounting for NPVf. The relationship between leaf pigment content and spectral indices varied depending on the scenario, due to the contribution of woody elements or brown pigments to the canopy reflectance. Therefore these statistical models aiming at estimating leaf pigment content integrated different mixtures of materials assumed to be found in experimental data (wood, senescent leaves).
Despite differences in these relationships, a strong correlation was obtained between all pigments and corresponding spectral indices for all scenarios accounting for NPVf in Site A. The main influence of NPVf was on the range of spectral indices, thus revealing different distributions of pigment content expected when applying Site A regressions to Site B. The poor suitability of some of these models (such as the models leading to negative anthocyanin content) highlights weaknesses in the integration of some of the factors influencing canopy reflectance into simulations.
Further work would be required in order to perform direct validation of these models for the estimation of leaf pigments based on imaging spectroscopy. However, it was expected that such straightforward approach based on spectral indices would lead to inaccurate estimation when applied on such complex ecosystems. The main purpose of our approach was to provide a basis for the integration of relative variations in pigment content on our study site, and more refined multivariate approaches may be needed to minimize absolute bias in simulated values, that we could not control or correct here.
We compared the estimated pigment contents following NPVf scenarios with those reported in the literature in the context of tropical forests in order to determine their realism for our purposes. Arellano et al. [96] reported Cab ranging from 19.2 to 105.4 µg.cm -2 and Cxc ranging from 1.0 to 44.2 µg.cm -2 in an Amazon rainforest of Ecuador. The ranges obtained for the different simulations then remained in a realistic range for Cab and Cxc. Canth showed discrepancies in some cases, with negative values. Ranges for Canth reported in the literature can reach higher values than those derived from our hybrid models [63,84], but such high values correspond to leaves which are unlikely to be found as main contributors to canopy reflectance in tropical forests. In tropical forests, anthocyanins mainly accumulate in juvenile leaves during the expansion phase and in vegetation found in the shaded parts of forest. Further investigations on the definition of optimal NPVf could be relevant for improving the matching between estimated and expected pigment contents in tropical forests.
Despite the ubiquity of spectral indices in remote sensing, there are some issues when upscaling from leaf to canopy level since other factors such as canopy architecture, solar geometry, background properties and LAI also affect the relationship between pigment content and vegetation index [97]. Intensive sampling in situ combined with Partial Least Square (PLS) method, which has proven a successful empirical approach to retrieve leaf traits from canopy spectral data [98][99][100], could be relevant for local calibration and improvement of pigment content estimation.

Influence of NPVf on simulation of canopy reflectance
The consideration of woody elements within forest architecture is one of the main issues in forest RTM [49]. In this respect, all DART NPVf scenarios were compared to imaging spectroscopy in order to investigate the influence of NPVf on simulation of canopy reflectance. Particularly important influence of NPVf was observed in the blue and red parts of the VIS domain, corresponding to the main domains of absorption of photosynthetic pigments with a major contribution of Cab. Pigments are major contributors of pixel's reflectance in the visible region [37,101]. Our results showed that the introduction of NPVf from WAD (NPVfWAD and NPVfWAD/2) resulted in an increase of reflectance in the blue and the red domains (Fig. 8). This observation is due to the mixture effect of green vegetation (strong absorption in the blue and red domains) with nonphotosynthetic vegetation (lower absorption in the visible) which leads to an increase of reflectance. These findings were in agreement with Malenovský et al. [54] who evidenced an increase of 2% of canopy reflectance in the red domain when woody elements were integrated into simulations. The introduction of NPVf based on brown pigments only resulted in the addition of absorbing constituents to the leaves in the simulation, and no mixture with elements of higher reflectance. Consistent with previous observations [54,102], the introduction of NPVf led to a decrease of simulated reflectance in the red edge and NIR, either using WAD or Cbrown. Overall, we found lower agreement between imaging spectroscopy and DART simulations in the NIR.
The reflectance in the red region was underestimated for NPVf0 and our optimal scenario (NPVfCbrown) compared to their counterpart in the imaging spectroscopy. Defining NPVf as a fraction of PAD (NPVfWAD) needs further attention. Combining the integration of woody elements through WAD and senescent materials through brown pigments in simulations may further improve the agreement between imaging spectroscopy and DART simulations.
Two assumptions made when considering NPVf as a uniform fraction of PAD should be discussed. First this leaf-to-wood area ratio is known to be spatially heterogeneous (across species, with canopy depth, etc.). Second the turbid medium assumption is expected to lead to discrepancies when introducing woody elements which are spatially hierarchically organized and interact with light as volumes instead of surfaces. The additional implicit assumption of homogenous mixing between leaf and wood is also not met. Improved representation of tree branching structure based on connected cylinders may contribute to improve the integration of woody elements [103,104]. A mixed representation using voxels for foliage and meshed volumes for woody parts is entirely compatible with the DART model, however the ALS point cloud alone does not provide sufficient resolution to achieve such level of detail of canopy geometry description. Leaf and wood are virtually inseparable in ALS point clouds as many returns are mixed points (with both leaves and wood contributing to a particular return) due to the laser footprint size (around 25cm diameter here). Tree geometry reconstruction would further be hampered by the limited pulse density, the geometric accuracy (limited by the laser footprint size), and the high occlusion rate occurring in dense forests. Finally, small woody elements such as the twigs have been shown to significantly affect the reflectance more than large woody parts such as branches or trunks [54], which suggests that the improved integration of the main wood structure in a mixed representation, which could be considered with TLS measurements, should also be refined with a procedure allocating small elements.

Influence of LOP variability on species separability and discrimination
The inadequate representation of LOP spatial variability within crowns, among crowns or species was reported as one of the main sources of error in RTM studies [49,105,106]. To solve these issues, we tested three scenarios accounting for LOP spatial variability : LOPspecies, LOPITC and LOPpixel.
Agreement between imaging spectroscopy and DART simulations was assessed by comparing the whole spectrum using statistical metrics, species separability using spectral dissimilarity metrics and potential of species discrimination using discriminant analysis. The lowest agreement was obtained by LOPspecies representation. This scenario led to overestimation of species separability. The LOPITC representation showed intermediate results while the most faithful scenario was LOPpixel.
The efforts made in the definition of LOPs turned out to be an important discriminating axis of species. But, we were constrained by the lack of detailed description of other sources of variability of the LOPs such as the foliar clumping, N, LADF or LMA homogenized on the whole 3D mockup. This potentially explain the lower discrimination of species in simulated images (LOPpixel and LOPITC).
The fine scale definition of leaf optical properties partly allows conservation of the spatial heterogeneity of reflectance, and appears as a promising framework exploring the Spectral Variation Hypothesis by means of remote sensing simulations.

Conclusion
DART RTM was used to simulate imaging spectroscopy in tropical forests. We took advantage of ALS to reconstruct a detailed forest 3D architecture. With the aim of defining a framework for producing realistic simulation in tropical forest, the influence of NPVf and LOP variability on simulated reflectance was explored. Analyses were performed using sunlit pixels of a set of ITCs. In order to estimate pigment content used as input for LOP characterization of each pixel, we used a hybrid approach combining 3D RTM and spectral indices.
NPVf were incorporated into simulations either considering WAD or Cbrown as proxies. Three scenarios including LOPspecies, LOPITC and LOPpixel enabled the assessment of the effect of LOP variability on simulated reflectance. Comparison between imaging spectroscopy and DART simulated images was performed using statistical metrics, spectral dissimilarity metric and supervised classification.
Our results show that simulations which did not account for NPVf resulted in moderate correlation (R 2 0.10-0.62) between measured and simulated canopy reflectance in the VIS domain, and consistently low correlation (R 2 0.08-0.18) in the red edge and NIR domain, associated with relatively high RMSE. Using Cbrown as a proxy for NPVf resulted in similar performances in the VIS domain, lower RMSE between measured and simulated canopy reflectance in the red edge and NIR domains, but still relatively low correlation. Since good agreement between imaging spectroscopy and DART simulation were obtained for NPVfCbrown, the underestimation of the reflectance in the blue and the red highlighted the need of combining the two approaches by defining an optimum NPVfWAD.
Our results showed that defining specific LOP for each pixel in a DART simulation resulted in more realistic simulations than when using the same LOP for all pixels in a tree crown, or when using the same LOP for all pixels in all tree crowns of a species. This was observed for radiometry, spectral dissimilarity and performances when applying supervised classification. We conclude that the strategy usually adopted in the literature, consisting in assigning one set of LOP per species results in strong overestimation of the capacity of supervised approaches to perform species discrimination. Spectral variation hypothesis needs to consider the large within species variability when inferring species diversity patterns from spectral diversity.
Overall, our results show that the incorporation of NPVf combined with LOP variability are relevant for simulating spatial heterogeneity in remotely sensed data in tropical forests. However, the improvement of leaf trait definition and characterization of structural parameters should be addressed in order to increase radiometric agreement between simulation and imaging spectroscopy data. Our simulation framework showed interesting potential for benchmarking remote sensing applications dedicated to biodiversity monitoring when applied to satellite missions currently operational as well as for future missions currently in preparation [107][108][109].