Simulating solar-induced chlorophyll fluorescence in a boreal forest stand reconstructed from terrestrial laser scanning measurements

Solar-induced chlorophyll fluorescence (SIF) has been shown to be a suitable remote sensing proxy of photosynthesis at multiple scales. However, the relationship between fluorescence and photosynthesis observed at the leaf level cannot be directly applied to the interpretation of retrieved SIF due to the impact of canopy structure. We carried out a SIF modelling study for a heterogeneous forest canopy considering the effect of canopy structure in the Discrete Anisotropic Radiative Transfer (DART) model. A 3D forest simulation scene consisting of realistic trees and understory, including multi-scale clumping at branch and canopy level, was constructed from terrestrial laser scanning data using the combined model TreeQSM and FaNNI for woody structure and leaf insertion, respectively. Next, using empirical data and a realistic range of leaf-level biochemical and physiological parameters, we conducted a local sensitivity analysis to demonstrate the potential of the approach for assessing the impact of structural, biochemical and physiological factors on top of canopy (TOC) SIF. The analysis gave insight into the factors that drive the intensity and spectral properties of TOC SIF in heterogeneous boreal forest canopies. DART simulated red TOC fluorescence was found to be less affected by biochemical factors such as chlorophyll and dry matter contents or the senescent factor than far-red fluorescence. In contrast, canopy structural factors such as overstory leaf area index (LAI), leaf angle distribution and fractional cover had a substantial and comparable impact across all SIF wavelengths, with the exception of understory LAI that affected predominantly far-red fluorescence. Finally, variations in the fluorescence quantum efficiency (Fqe) of photosystem II affected all TOC SIF wavelengths. Our results also revealed that not only canopy structural factors but also understory fluorescence should be considered in the interpretation of tower, airborne and satellite SIF datasets, especially when acquired in the (near-) nadir viewing direction and for forests with open canopies. We suggest that the modelling strategy introduced in this study, coupled with the increasing availability of TLS and other 3D data sources, can be applied to resolve the interplay between physiological, biochemical and structural factors affecting SIF across ecosystems and independently of canopy complexity, paving the way for future SIFbased 3D photosynthesis models. https://doi.org/10.1016/j.rse.2019.111274 Received 28 July 2018; Received in revised form 3 June 2019; Accepted 19 June 2019 ⁎ Corresponding author at: LREIS, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China. E-mail addresses: lww@lreis.ac.cn (W. Liu), jon.atherton@helsinki.fi (J. Atherton), jean-philippe.gastellu@iut-tlse3.fr (J.-P. Gastellu-Etchegorry), pasi.raumonen@tut.fi (P. Raumonen), markku.akerblom@tut.fi (M. Åkerblom), raisa.makipaa@luke.fi (R. Mäkipää), joan.porcar@helsinki.fi (A. Porcar-Castell). Remote Sensing of Environment xxx (xxxx) xxxx 0034-4257/ © 2019 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/). Please cite this article as: Weiwei Liu, et al., Remote Sensing of Environment, https://doi.org/10.1016/j.rse.2019.111274

Solar-induced chlorophyll fluorescence (SIF) has been shown to be a suitable remote sensing proxy of photosynthesis at multiple scales. However, the relationship between fluorescence and photosynthesis observed at the leaf level cannot be directly applied to the interpretation of retrieved SIF due to the impact of canopy structure. We carried out a SIF modelling study for a heterogeneous forest canopy considering the effect of canopy structure in the Discrete Anisotropic Radiative Transfer (DART) model. A 3D forest simulation scene consisting of realistic trees and understory, including multi-scale clumping at branch and canopy level, was constructed from terrestrial laser scanning data using the combined model TreeQSM and FaNNI for woody structure and leaf insertion, respectively. Next, using empirical data and a realistic range of leaf-level biochemical and physiological parameters, we conducted a local sensitivity analysis to demonstrate the potential of the approach for assessing the impact of structural, biochemical and physiological factors on top of canopy (TOC) SIF. The analysis gave insight into the factors that drive the intensity and spectral properties of TOC SIF in heterogeneous boreal forest canopies. DART simulated red TOC fluorescence was found to be less affected by biochemical factors such as chlorophyll and dry matter contents or the senescent factor than far-red fluorescence. In contrast, canopy structural factors such as overstory leaf area index (LAI), leaf angle distribution and fractional cover had a substantial and comparable impact across all SIF wavelengths, with the exception of understory LAI that affected predominantly far-red fluorescence. Finally, variations in the fluorescence quantum efficiency (Fqe) of photosystem II affected all TOC SIF wavelengths. Our results also revealed that not only canopy structural factors but also understory fluorescence should be considered in the interpretation of tower, airborne and satellite SIF datasets, especially when acquired in the (near-) nadir viewing direction and for forests with open canopies. We suggest that the modelling strategy introduced in this study, coupled with the increasing availability of TLS and other 3D data sources, can be applied to resolve the interplay between physiological, biochemical and structural factors affecting SIF across ecosystems and independently of canopy complexity, paving the way for future SIFbased 3D photosynthesis models.

Introduction
Solar-induced chlorophyll fluorescence (SIF) is electromagnetic radiation emitted by plants during daylight in the red and near-infrared wavelengths. Measurable from space and airborne platforms, SIF originates from within the photosynthetic apparatus of all higher plants and is directly linked to carbon assimilation and plant physiological status (Porcar-Castell et al., 2014). Unlike traditional greenness vegetation indices calculated from reflectance, such as the Normalized Different Vegetation Index (NDVI), SIF responds very rapidly to vegetation stress, and therefore has the potential to track time-varying plant functional dynamics both at seasonal and (sub) diurnal scales (Amoros-Lopez et al., 2008;Fournier et al., 2012;Yang et al., 2015;Alonso et al., 2017).
The magnitude of SIF is very small compared to the radiation reflected by plant canopies being less than 5% of reflected sunlight in the near-infrared (Meroni et al., 2009). Additionally, under natural illumination, SIF cannot be measured directly due to the spectral overlap with reflected radiation. Instead, SIF is estimated within Fraunhofer or telluric absorption bands, such as the O 2 -A band around 761 nm and O 2 -B band around 687 nm using the Fraunhofer Line Depth (FLD) technique . At the leaf scale, spectral fluorescence is not only influenced by the dynamics of photosynthesis, but also by chlorophyll a + b (Cab) content and leaf structure (Van Wittenberghe et al., 2015;Atherton et al., 2017;Magney et al., 2019a). In addition, at the canopy scale, SIF is influenced by 1) sun and view geometry, 2) canopy structure, 3) instrumental effects, and 4) the selection of SIF retrieval algorithm (Guanter et al., 2010;Damm et al., 2015a;Liu et al., 2016).
In recent years, SIF has been used to characterize the spatiotemporal dynamics of photosynthetic gross primary production (GPP), a key variable for carbon cycle and climate change studies (Frankenberg et al., 2011;Guanter et al., 2014;Zhang et al., 2014;Sun et al., 2017Sun et al., , 2018Magney et al., 2019b). The relationship between SIF and GPP in these studies is likely to be at least partially affected by the spatial and temporal variation in vegetation structure across biomes. In a modelling study, Verrelst et al. (2015) estimated that approximately 80% variation of TOC SIF can be attributed to the leaf optical properties and canopy structure. In contrast, recent studies suggest that canopy structure effects might to some extent cancel out at the scale of a satellite pixel when relating SIF to GPP (Sun et al., 2017;. Hence, more research is required to disentangle the complex relationship between SIF and GPP. In particular, schemes for coupling leaf fluorescence models with canopy radiative transfer (RT) contribute towards an improved interpretation of SIF observation in terms of plant physical, biochemical and physiological traits.
Early attempts at using physically based models to simulate SIF were carried out by the European Space Agency (ESA), in which a leaf biophysical model FluorMODleaf and leaf-canopy fluorescence model FluorSAIL were developed (Miller et al., 2004) based on the Scattering of Arbitrarily Inclined Leaves (SAIL) RT model (Verhoef, 1984). Later, a more comprehensive model system, SCOPE (Soil-Canopy Observation, Photosynthesis and Energy Balance) was developed as part of the early Fluorescence Explorer (FLEX) mission preparatory studies. SCOPE coupled the turbid medium SAIL with the Soil-Vegetation-Atmosphere (SVAT) models in order to simulate photosynthesis, thermal radiation and energy balance, including SIF ( Van der Tol et al., 2009). However, since SAIL is a one-dimensional (1D) RT model, designed for homogenous vegetation canopies, it is not suitable for assessing the impact of horizontal and vertical heterogeneity in structurally complex canopies such as forest (Porcar-Castell et al., 2014;Migliavacca et al., 2017). Recently, a new multilayer version of the SCOPE model (mSCOPE) was developed to partially address this shortcoming by simulating vertically heterogeneous canopies (Yang et al., 2017).
Unlike 1D models, 3D RT models have been developed to address the need for SIF simulations of structurally complex canopy architectures. A coupled FluoMODleaf-FluorSAIL model was improved to simulate the SIF in heterogeneous canopies with the introduction of the first-order approximation forest model (FLIM) (Zarco-Tejada et al., 2013). Yet, crown scale specific effects due to the interplay of sunlit and shaded crown fractions, caused by crown element overlap, shadowing and light scattering, were not considered in FLIM. More recently, Zhao et al. (2016) introduced FluorWPS, which is a Monte Carlo ray-tracing based 3D radiative transfer model for row-planted canopies, such as soybean and cotton. Another radiative transfer model, FluorFLIGHT  couples the 3D ray-tracing forest radiative transfer model FLIGHT and the leaf fluorescence model FLUSPECT. FLIGHT describes the tree canopy as consisting of homogeneous tree crowns and gaps between them (North, 1996). The results of FluorFLIGHT showed that the TOC SIF signal is greatly influenced by canopy structure for complex canopy types with heterogeneity in both horizontal and vertical dimensions. Although these two models represent a significant step forward in complexity, they approximate plant crowns as geometrical object filled homogeneously with leaves without considering leaf clumping effects. Contrary to this, the 3D Discrete Anisotropic Radiative Transfer (DART) model can be used to provide a more realistic representation of plant architecture by explicitly describing the foliage distribution, including leaf clumping at branch and crown levels, along with detailed geometry of stems and branches (Gastellu-Etchegorry et al., 2017). DART also simulates spectrally resolved radiative budget of vegetation scenes, which can be used to directly estimate the light regime and instantaneous incoming photosynthetically active radiation (PAR) of each leaf. Subsequently, it is possible to separate foliage into sun and shade adapted parts and assign them specific optical properties and fluorescence yields as a function of light regime or incoming PAR (Gastellu-Etchegorry et al., 2018). The parameterization of leaf geometry, distribution and clumping at branch and crown levels can, however, be very laborious and consequently forced to be simplified, especially for a large canopy.
In this study we present a methodological scheme to assimilate high resolution forest canopy 3D structural data, acquired with a terrestrial laser scanner in a Silver birch stand, into the DART model. The scheme is combined with empirical data and used to simulate TOC SIF in the heterogeneous forest stand with consideration to multiscale leaf clumping and understory vegetation. The potential of the modelling scheme is next demonstrated by conducting a local sensitivity analysis of TOC SIF to key structural, biochemical and physiological factors.

Site description
The simulated scene was based on a TLS surveyed circular forest plot (25 m in radius) dominated by Silver birch (Betula pendula Roth) and located in the vicinity of the Station for Measuring Ecosystem-Atmosphere Relations II (SMEAR II), in southern Finland. The coordinates of the center-point of the circular site were 24.31478°E, 61.84335°N. The plot included 185 Silver birch trees and 5 Norway spruce trees. At the peak growing season, the understory layer was composed of grasses, shrubs and bryophytes.

Terrestrial laser scanner data
Terrestrial laser scanning (TLS) measurement survey was carried out in March of 2017 before leaf burst. A Leica ScanStation P40 (Leica Geosystems, Heerbrugg, Switzerland) was used to perform the point cloud scanning. Five scans were carried out in different positions within the plot to cover most of the trees and to minimize the influence of occlusion of trees: one scan was performed in the center of the plot, and four scans at the edge of the plot with a cross sampling strategy. Reflectors were attached across each tree perimeter within the plot at breast height (1.3 m), and used to facilitate the registration of the point W. Liu, et al. Remote Sensing of Environment xxx (xxxx) xxxx cloud scanned from five viewpoints. The point cloud dataset was collected initially with 3.1 mm spacing at 10 m distance and then downscaled to 5 mm spacing at 10 m using Leica Cyclone 9.1.4 (Leica, 2017) to reduce the data size. Wind speed was very small (below 3 m/s) during the measurements and the influence caused by the moving of crown elements was not considered in this paper.

Digital hemispherical photographs
We used the HemiView Forest Canopy Image Analysis System (Delta-T Devices Ltd, Cambridge, UK) to obtain the upward digital hemispherical photographs (DHP) of the birch forest. This canopy image analysis system was composed of a Canon EOS 70D camera equipped with a fisheye lens (Sigma 4.5 mm f/2.8 EX DC), and a special gimbal stabilizer with a compass that keeps the camera pointing upwards and mark the North and South in the resulting image. The canopy image analysis system was mounted on a tripod with a height of 1 m. A total of 9 digital hemispherical photographs were collected in the study plot on the same day as the terrestrial laser scanner data and using a cross-shaped sampling strategy. The measured DHP pictures were then processed using the CAN-EYE V6 (Weiss and Baret, 2010) to calculate the gap fraction and LAI.

Spectral data
Leaf reflectance and transmittance data were collected during the summer of 2014 using an ASD hand held spectrometer (Malven Panlytical Inc., Westborough, USA) and ASD integrating sphere (PANlytical Inc., Westborough, USA). The reflectance spectra of birch bark and understory were measured at the plot using the same spectrometer during summer of 2017. All spectral data were measured with the spectral resolution of 3 nm in the spectral range of 400 to 1000 nm. Leaf reflectance and transmittance data were used for the estimation of leaf structure parameter (N). The reflectance spectra of birch bark and understory were used for reflectance and SIF modelling in DART. Spectral data are shown in Fig. S1  Owing to the limited number of TLS scans, constructing each single tree within the experimental plot was not feasible due to occlusion effects. Therefore, a point cloud dataset of a single tree that was least influenced by occlusion was used to reconstruct a single quantitative tree structure model (QSM). We used TreeQSM (Raumonen et al., 2013;Calders et al., 2015) and a non-intersecting leaf insertion algorithm (FaNNI) (Åkerblom et al., 2018) to reconstruct the woody and leafy elements of the birch tree, respectively.
Constructing the woody elements: to construct the woody elements of a birch tree, the laser scanner data were processed using TreeQSM. In TreeQSM, the TLS data is segmented into individual woody elements (stem and branches) after applying a patchwork 'cover sets' approach which grows a global surface by connecting local patches (Raumonen et al., 2013). In the latter steps of the method, the segmented surface is fitted with cylinders to obtain the final compact (relative to the TLS data) tree model. Leaf shape parameterization: we used the same leaf size and shape used for the fourth phase of the Radiation transfer Model Intercomparison (RAMI) at Järvselja birch stand, Estonia. The parameterized birch leaf was made up of eight triangles to depict the main outline of birch leaf (Fig. S2 in the Supplementary material). A petiole of 5 cm was generated to connect it to the branch. Foliage insertion: a non-intersecting leaf insertion algorithm (FaNNI), not allowing leaves to intersect with each other and other elements, was used to perform the leaf insertion process based on the generated woody elements. Leaves were distributed around the ends of the branches generated during the woody element construction. In other words, the leaf distribution within the crown was indirectly informed by the laser scanner data. The phyllotaxy, the arrangement of leaves on stem and branch, was assumed to follow the Fibonacci sequence: the angle of adjacent leaves in the plane perpendicular to the branches was 137.5°( Fig. S3 in the Supplementary material). Fibonacci sequence has been previously used also for other structural levels, for example, within a shoot (Smolander and Stenberg, 2003;Disney et al., 2006). Leaf angle distribution: A planophile leaf angle distribution with a mean leaf angle of 31.04 was adopted here as a stand-average value for birch (Betula pendula) which was measured at Bergius Botanical Garden of Sweden using levelled-digital camera approach for the same tree species (Pisek et al., 2013). These studies indicated that a planophile or plagiophile leaf angle distribution (LAD) is more appropriate for the deciduous broadleaf tree than the commonly used spherical distribution. The exact number of leaves and their insertion points were only determined after all trees in the stand were generated as described below.
The reconstructed 3D tree is shown in Fig. 1. More detailed biometric information can be found in Fig. S4 in the Supplementary material.

Creating a structurally complex 3D forest scene
We determined the precise coordinates and dimensions of each tree in the scene from point cloud data using 3D forest software (Trochta et al., 2017). For each tree, we calculated three scale factors by dividing each tree dimension (tree height, crown width and depth) with the corresponding value for the model tree. Hence the virtual forest scene ( Fig. 2) was populated with scaled replicates of the model tree positioned at the appropriate coordinates. It should also be noted that leaf size was kept constant across the scene.
The number of leaves in a tree was needed to determine their insertion points along the branches. We developed an optimization algorithm to obtain the threshold (leaf number) to determine the leaf insertion process. We simulated a series of forest scenes with decreasing distances between adjacent leaves (i.e., increasing number of leaves and LAI) with a LAI step of 0.1. The number of leaves was then optimized by comparing the gap faction of the simulated stand to that retrieved from hemispherical photographs. Subsequently, LAI of the reconstructed 3D forest scene could be calculated from the number of leaves and single leaf area.
Once completed, the digitized stand was used as the input scene for the canopy radiative transfer model DART, which was coupled to the leaf fluorescence model FLUSPECT. The complete scene construction and SIF simulation framework is outlined as a flowchart in Fig. 3.

SIF modelling and local sensitivity analysis
In this section we describe the parameterization of FLUSPECT inside DART, as well as the range of values selected for the local sensitivity analysis.

Modelling of leaf fluorescence using FLUSPECT
DART uses the FLUSPECT model (Vilfan et al., 2016), which is based on the widely used PROSPECT model (Jacquemoud et al., 2009), as the fluorescence source within our modelling framework. FLUSPECT calculates the forward (abaxial) and backward (adaxial) SIF spectrum from 640 to 850 nm of photosystem I (PSI) and photosystem II (PSII). The magnitude of the SIF signal of PSI and PSII modeled by FLUSPECT is controlled by three main factors: (1) fluorescence quantum efficiency (Fqe) of PSI and PSII, (2) the amount of leaf absorbed photosynthetically active radiation (APAR), and (3) leaf structural and biochemical variables controlling light scattering and absorption inside the leaf.
The Fqe of PSI and PSII is the emission efficiency of chlorophyll a fluorescence at the photosystem level. It is defined as the probability of an absorbed photon to be re-emitted as fluorescence by PSI or PSII, respectively. As a first approximation, and since the fluorescence quantum yield of PSI is known to remain stable in response to shortterm illumination (Genty et al., 1990;Pfündel, 1998), Fqe of PSI was assumed to be constant. Accordingly, it was assumed that variations in the overall Fqe originate from variations in the fluorescence efficiency of PSII alone.
We need to emphasize that the spatial and temporal dynamics of PSI fluorescence remain very unclear, and therefore the assumption of constant PSI fluorescence used here should be taken only as a first approximation. In fact, this lack of understanding has motivated changes in the latest FLUSPECT CX parameterization (Vilfan et al., 2018), where both PSII and PSI fluorescence contributions are merged into a single component. For comparison, an example of the maximum range of variation in leaf level fluorescence in response to a tenfold change in Fqe, when using the two different FLUSPECT formulations, is presented. As illustrated in Fig. 4, both the shape of the spectra and its variation in response to Fqe differ substantially between versions. Although the new version of FLUSPECT appears to better reproduce the observed shape of fluorescence spectra of birch leaves (Rajewicz et al., 2019), we decided to use the previous version, as we felt that presenting a scheme that pools together PSII and PSI fluorescence would not contribute to promoting understanding of the possible functionality and dynamics of PSI fluorescence.   , 1996;Gastellu-Etchegorry et al., 2017) was used to perform upscaling of the SIF spectrum from leaves to canopy in the parameterized forest scene described above. DART is a physically based 3D radiative transfer model, which has been designed to simulate radiation propagation from visible to thermal infrared in complex and heterogeneous 3D landscape scenes, such as forest scenes (Malenovský et al., 2008).
A square scene was constructed with the constructed trees ( Fig. 2) which were simulated by geometrically explicit facets. Assuming that the simulated forest plot is a representative sample of forest, we used an option of repetitive scene to eliminate boundary effects, which means that the radiative flux escaping the scene from one of its vertical faces entered the scene at the corresponding position of the scene opposite side. Additionally, an understory layer (baseline LAI equal to 1.0) was created to simulate SIF emitted by the understory. The understory was simulated as a 30 cm high layer of triangular facets randomly distributed in space according to a spherical LAD. For simplicity, understory foliar elements were given the same optical and fluorescence parameters as the tree foliar elements. The optical properties of branches and twigs were given the same reflectance value as measured for bark. Top of atmosphere (TOA) sun irradiance is defined from the DART database. The bottom of atmosphere (BOA) direct solar and diffuse irradiances were simulated with "USSTD 76" gas model and the "USSTD Rural 23km" aerosol model (Fig. 4b). We used a solar position of 30°s un zenith and 225°sun azimuth angle, which corresponded to a summer noon at the experimental site. Table 2 shows the stand-average and extreme value used for DART model parametrization.
The fluorescence radiative transfer in DART was conducted with the flux tracking method, which propagates the forward-/backward- emitted (abaxial/adaxial) fluorescence simulated by the FLUSPECT leaf model into a finite number of discrete directions using an iterative scattering approach. TOC SIF was defined by the fluorescence radiance that escaped from the canopy in nadir viewing direction. It was modeled simultaneously with TOC reflectance from 400 nm to 850 nm using the constructed 3D forest scene and corresponding biochemical and physiological variables between 640 and 850 nm at a spectral resolution of 1 nm. In addition, radiation budget products (e.g. absorbed photosynthetically active radiation -APAR) were also simulated for each pixel in the nadir image and element type (ground, woody elements and leaves) as the sum of the energy absorbed by the elements at the spatial location corresponding to the pixel, regardless of them being visible to the sensor. The radiative budget of irradiance absorbed by leaves resulted in the actual foliage APAR, which was used to calculate a normalized TOC SIF. More detailed information about the model setup and the normalization of TOC SIF is available in the Supplementary material (see Supplement Sect. 4 and 5).

Sensitivity analysis
We explored the relative contributions of structural, biochemical and physiological factors to TOC SIF by means of a local sensitivity analysis (i.e. changing one variable for each simulation). We tested both leaf-biochemical and canopy-structural factors, and where possible our parameter values (ranges) were rooted in field observations.
Empirically estimated levels of Cab, Car, and dry matter obtained in the same experimental plot (Atherton et al., 2017), were used as default stand-average variables in the model. For the sensitivity analysis, the variation range for foliar pigment contents was selected on the basis of a recent global meta-analysis of plant pigments (Esteban et al., 2015). The default stand-average variables and variation ranges used in the sensitivity analysis are summarized in Table 1. The leaf structural parameter (N) was estimated by fitting the FLUSPECT modeled leaf reflectance and transmittance spectra, generated with the default stand-average biochemical variables (see Table 1), to measured leaf reflectance and transmittance spectra (see Fig. S1 in the Supplementary material). Based this estimation, the leaf structure parameter was kept constant (N = 1.78) throughout this study.
The fluorescence yield efficiency (Fqe) was estimated from PAM fluorescence measurements as 0.0154, 0.0201 and 0.0053 for PSII of sunlit leaves, PSII of shaded leaves and PSI of all leaves, respectively (see Supplement Sect. 3). For the sensitivity analysis, we used a PSII efficiency variation range from 0.0015 to 0.0154 for sunlit leaves and from 0.0020 to 0.0201 for shaded leaves, in an attempt to replicate the full natural range of physiological variation observed in steady state fluorescence quantum yield (F s ) where the lower end would reflect a situation with strongly downregulated foliage (Porcar-Castell et al., 2008).
The LAI value was modified using two approaches. First, when the LAI was changed from the stand-average value (LAI equal to 2.7), calculated in Section 2.2.2, to the minimum value (LAI equal to 0.75), it was decreased by decreasing leaf dimensions, without modifying the woody elements. This situation approximates, in the inverted temporal direction, the seasonal increase in LAI when leaves are unfolding and expanding in spring. Second, when the LAI was changed from standaverage to the maximum value (LAI equal to 6.0), it was increased by increasing the number of leaves while keeping leaf size constant. The fraction of tree cover (fCover) was modified from a sparse forest to a dense forest by adjusting the scene horizontal dimensions (i.e. tree coordinates) without modifying the number of trees, their size, nor their architecture. This approach would simulate changes in tree density, such as upon a commercial forest thinning.
Finally, two types of understory were created for the sensitivity analysis: (1) understory with no SIF emission and LAI equal to 1.0, (2) understory with SIF emission and LAI equal to 0.5, 1.0, 1.5, and 2.0.

Architectural characteristic of the 3D forest scene
Directional gap fractions derived from the digital hemispherical cameras and DART simulations for 6 zenith angle ranges are shown in Fig. 5. From 10°to 60°, the two calculations of gap fraction were consistent and had relatively small variation ranges, denoting that the constructed forest scene matched with the real forest scene in structural characteristics. The relative difference (error) between simulated and measured gap fractions was slightly larger at smaller zenith angles from 0°to 10°, as were the uncertainty estimates in both measurements and simulations.

Spectral characteristics of the 3D forest sub-scenes
The TOC reflectance and SIF images were produced using the default stand-average variables for structural, biochemical and physiological variables presented in Tables 1 and 2. We subjectively selected four representative areas, here called sub-scenes, from the false color synthesis image for further analysis: sunlit crown, shaded crown, sunlit understory and shaded understory (Fig. 6).

Reflectance spectra of the sub-scenes
As expected, the modeled reflectance spectra of the four sub-scenes (sunlit crown, shaded crown, sunlit understory and shaded understory sub-scene) were found to be different in shape and variation ranges (Fig. 7). The sunlit crown sub-scene (Fig. 7a) had the largest internal variation range due to the spectral heterogeneity originating from sunlit and shaded leaves. Note the peak in the apparent reflectance at 760 nm induced by fluorescence emission in all sub-scenes.

SIF spectra of the sub-scenes
The four forest sub-scenes also presented different SIF spectra (Fig. 8). Sunlit crown sub-scene had the maximal internal variation range compared to the other three sub-scenes (Fig. 8a). Interestingly, even though the understory had a smaller LAI (LAI = 1) compared to the overstory (LAI = 2.7), understory sub-scenes still emitted a considerable SIF signal (Fig. 8b, d).

Impact of structural factors on TOC SIF
TOC SIF spectra were simulated for several ranges of the structural factors: overstory LAI, LAD, fCover and understory LAI with SIF emission (Fig. 9). In general, variations in structural factors strongly affected the magnitude of TOC SIF across all wavelengths with only moderate changes in the shape of the spectra (Fig. 10). The red/far-red ratio tended to decrease with increasing LAI in a similar fashion for overstory and understory vegetation. In contrast, fCover and LAD had only a minor positive effect on the red/far-red ratio (Fig. 10c).
TOC SIF at the two emission peaks revealed that if overstory LAI increased from 0.75 to 6.0, the TOC SIF increased accordingly without saturation (Fig. 10). Spherical and vertical leaf angle distributions lead to a SIF decrease of about 25% and 60%, respectively, when compared with the planophile LAD (Fig. 9c). In relation to canopy coverage, a decrease in fCover (e.g. response to forest thinning or forest disturbance) resulted in almost linear decrease in TOC SIF (Figs. 9c, 10). Understory SIF had a larger effect on the TOC SIF in the far-red SIF region (740 nm) compared with the more modest increase in the red SIF region (685 nm), consistent with the re-absorption of red fluorescence by the canopy above. Importantly, leaves in understory made a considerable contribution to TOC SIF when compared to leaves in overstory (Fig. 10). For instance, in the far-red region, increasing understory LAI from 0.0 to 2.0 led to a TOC SIF increase of approximately 0.9 mW/m 2 / sr/nm, while the same TOC SIF increase could only be achieved through a larger increase in overstory LAI (e.g., overstory LAI increase from 0.75 to 6.0). This observation is clearly affected by the constrained canopy architecture during LAI change from 0.75 to 6.0, however it emphasizes the strong influence that understory SIF can have in open canopies.

Impact of biochemical and physiological factors on TOC SIF
TOC SIF spectra were also simulated for five biochemical variables: chlorophyll content, carotenoid content, water content, dry matter and senescent factor, and varying fluorescence quantum yields (Fig. 11). All biochemical variations, except the water content, had an effect on TOC SIF. Variations in biochemical and physiological factors not only affected the magnitude of TOC SIF but also its spectral shape (Fig. 12c).
Far-red SIF showed a remarkable positive correlation with increasing chlorophyll content (Fig. 11a) and a negative correlation with dry matter content (Fig. 11b) and senescent factor (Fig. 11e). In contrast, red SIF instead was only slightly and negatively affected by these three factors (Figs. 11, 12a). Carotenoid contents negatively affected TOC SIF across wavelengths (Fig. 11d). These results suggest that TOC far-red SIF (740 nm) tends to be more sensitive to the biochemical parameters than red SIF (685 nm). As expected, the Fqe of photosystem II similarly affected all wavelengths (Figs. 11,12) with SIF approaching zero with decreasing Fqe in the red bands but not in the far-red, due to  our preliminary assumption of constant and invariable PSI fluorescence (Fig. 4a).

Impact of chlorophyll content on TOC SIF and TOC SIF normalized by APAR
To expand on the mechanism by which changes in foliar chlorophyll content affect the TOC SIF we characterized the impact of chlorophyll content on both canopy level APAR (estimated by integrating the irradiance absorbed by leaves in the 400-700 nm range and across the vertical profile of canopies), and the relative fluorescence yield SIF/ APAR (Fig. 13). As expected, increasing chlorophyll content increased canopy absorption especially in the green and red wavelengths but not in the blue wavelengths which tended to be more efficiently absorbed by chlorophyll even at small chlorophyll content concentrations (Fig. 13c). The effect of chlorophyll content on the apparent TOC SIF/ APAR was clearly different from that of TOC SIF, presenting a reversed pattern of response to chlorophyll content in red and far-red SIF (Fig. 13b). In particular, the apparent yield of red SIF decreased with increasing chlorophyll content under the action of increasing re-absorption. In contrast, the far-red SIF yield, which is much less affected by re-absorption, was found to increase with increasing chlorophyll content. This increase can be explained as the gradual response of farred SIF to the increasing chlorophyll content contribution to total APAR relative to the background PAR absorption by carotenoids, dry matter and senescent material. In summary, chlorophyll content has a minor effect on TOC red SIF, because the positive effect of chlorophyll content on APAR and the negative effect on reabsorption partly cancel out. In contrast, chlorophyll content has a large and positive effect on TOC far- Fig. 6. a) False color composite image with a spatial resolution of 0.125 m simulated by DART using the default stand-average variables. Four sample areas refer to sunlit crown sub-scene (A), shaded crown sub-scene (B), sunlit understory sub-scene (C), and shaded understory sub-scene (D). b) The corresponding SIF image at 685 nm. c) The corresponding SIF image at 740 nm. The false color composite was simulated with bare soil background (i.e., no understory) to distinguish different sub-scenes. Understory with SIF emission and LAI = 1 was used to simulate SIF at 685 nm and 740 nm. red SIF through APAR (Fig. 13d).

Discussion
This study used TLS data to parameterize the fine 3D structure of a natural birch stand that was assimilated into the DART model. DART, coupled with the leaf fluorescence FLUSPECT model (Gastellu-Etchegorry et al., 2017), was parametrized to yield a quantitative scheme that could be applied to investigate the interaction between SIF, structural and functional properties found within a forest ecosystem. The potential of the scheme was demonstrated with a local sensitivity analysis which revealed some new insights into the factors that control TOC SIF in a complex forest ecosystem.

The realisation of a 3D forest scene constructed from TLS data for SIF modelling
We used TLS data and two open source processing algorithms TreeQSM (Raumonen et al., 2013;Calders et al., 2015) and FaNNI (Åkerblom et al., 2018) to model woody elements and foliage, respectively, and create a "virtual forest" (Fig. 2). We populated our forest scene with scaled realistic tree models based on a single tree construction (Fig. 1), then leaf area was tuned to match total canopy gaps. The realism of our virtual forest scene was demonstrated by the angular dependence of gap fraction in Fig. 5. The gap fraction coincided for the real and virtual stand. We also compared the DART simulated canopy spectral reflectance with a corresponding airborne imaging spectroscopy data (Markiet et al., 2017;Supplement Sect. 6) and found them to be within 0.058 reflectance units across the spectrum. We consider this as a reasonable match, owing to a number of unknown stand parameters (e.g., understory spectrum, leaf angle distribution, leaf scattering phase function, cloning of trees inside the scene, etc.). It is important to note that although we used a single tree model (due to the limited number of TLS acquisitions and occlusion effects), a full construction of every single tree is feasible in future work given a higher density of TLS sampling points (Calders et al., 2018).
Using a realistic characterization of our forest stand of interest, we investigated the spatial variability of SIF and reflectance within the simulated 3D scene containing both tree crown and understory (Figs. 7,8). As expected, the sunlit crown emitted substantially more fluorescence than the shaded crown and presented higher variation (Fig. 8). These phenomena were potentially caused by the heterogeneity of canopy structure and irradiance conditions within the tree crown (Damm et al., 2015b;Hernández-Clemente et al., 2017 andCamino et al., 2018), which relates to the spatial distribution of woody elements and foliage (Kuusk, 2018). In future studies, 3D quantitative models capable of reproducing fine variations in SIF within a plant crown could be used in connection with high spatial resolution data acquired with a SIF imaging system (Albert et al., 2019;Rossini et al., 2015). Such a combination could facilitate the spatially explicit investigation of plant stress in heterogeneous tree canopies (Camino et al., 2018).
The preponderance of natural gaps in the forest scene revealed an important and frequently overlooked aspect of canopy structure that has direct relevance to larger scale SIF observations; the importance of the vegetation understory layer. Several studies have found that the TOC reflectance signal is influenced, in extreme cases even dominated, by the understory vegetation and that this influence changes during the growing season (Eriksson et al., 2006;Kuusk, 2001;Kuusk et al., 2004;Rautiainen et al., 2007Rautiainen et al., , 2011Disney et al., 2011). For SIF, we found a similar level of importance for understory across our simulated scene. With a fCover of 70%, tree LAI of 2.7 and constant crown biochemical and physical variables, TOC far-red SIF value doubled when simulating an understory with LAI equal to 2 compared to the situation with no understory SIF (Figs. 9d and 10a, b). According to Pisek et al. (2015) and Liu et al. (2017), the understory LAI of a birch site may increase from zero to 2-3 from June to August. Importantly, a spatially-temporal responsive understory (in terms of TOC SIF) could have important consequences for the use of SIF in the estimation of canopy traits. Specifically, if the relative contribution of understory SIF to TOC SIF is not equivalent to the relative contribution of understory GPP to canopy GPP (Sakai et al., 2006;Kolari et al., 2006;Lin et al., 2018), the estimation of GPP from space could be significantly biased in open multistory ecosystems, which is an issue that calls for further investigation.

Leaf and canopy scale local sensitivity to natural variability in model parameters
We conducted a local sensitivity analysis as an investigation of multiscale SIF controls in our virtual scene. In absolute numbers, red SIF responded similarly to variations in Fqe than far-red SIF (Figs. 11f,12), but their relative range was different. The far-red TOC SIF did not converge to zero at Fqe equal to 0 (Fig. 12), which results from our assumption of a fixed PSI contribution. As demonstrated in Fig. 4, we would not achieve the same result if using the latest version of the FLUSPECT CX leaf model (Vilfan et al., 2018) that pools together the fluorescence components from both photosystems in a single Fqe parameter. Regardless of the modelling challenge, non-zero Fig. 9. Sensitivity of simulated forest scene TOC SIF spectra to changes in canopy and understory structural variables: a) overstory LAI, b) fCover, c) LAD, and d) understory LAI. Lines corresponding to increased variable values are in red, and lines for reduced variables are in blue. O_LAI refers to overstory LAI, U_LAI refers to understory LAI. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 10. The variability of TOC SIF caused by canopy structural variables and biophysical variables at two SIF peaks: a) TOC SIF in the red region (685 nm), b) TOC SIF in the far-red region (740 nm) and c) ratio of fed SIF to far-red SIF. O_LAI refers to overstory LAI, U_LAI refers to understory LAI. W. Liu, et al. Remote Sensing of Environment xxx (xxxx) xxxx fluorescence levels have been observed in boreal and alpine evergreens during periods when photosynthesis is fully inhibited by temperature (Ottander et al., 1995;Ensminger et al., 2004;Zarter et al., 2006;Porcar-Castell et al., 2008;Porcar-Castell, 2011). This suggests that a background fluorescence level is expected when the photosynthetic machinery is completely shut down. We, therefore, advocate for the use of models that explicitly separate PSII and PSI components, which would allow us to test and possibly better understand the sensitivity of model simulations to PSI fluorescence.
The local sensitivity analysis revealed that red SIF was less sensitive W. Liu, et al. Remote Sensing of Environment xxx (xxxx) xxxx to variation in chlorophyll content than far-red SIF (Fig. 12a, b). This result can be explained by two phenomena that affect different regions of the fluorescence spectra differently: i) increase in APAR that increases fluorescence across all wavelengths in response to increasing chlorophyll content, and ii) reabsorption that decreases fluorescence around the red wavelengths in response to increasing chlorophyll content. The combined result of these two processes was characterized by Gitelson et al. (1998) through experiment with ethanol chlorophyll suspensions of different concentration. At low but increasing chlorophyll concentrations, where SIF reabsorption is small and APAR increases linearly, both red and far-red fluorescence will increase with chlorophyll content. In case of increasing high chlorophyll concentrations, reabsorption of red fluorescence is larger relative to the increase in SIF emission by enlarged absorption of photosynthetically active radiation. This results in a net decrease in red fluorescence, while preserving an increase of far-red fluorescence. Importantly, the concentration at which the compensation point between APAR increase and SIF reabsorption occurs depends on the distribution and packing of the chlorophyll molecules, and ultimately on the spatial scale of the measurement (Romero et al., 2018). Accordingly, while our observed patterns are consistent with published results from chlorophyll suspensions (Gitelson et al., 1998;Buschmann, 2007) and leaf-level simulations , they suggest that at the canopy level, the point at which SIF reabsorption starts to compensate for increase in APAR takes place at very low chlorophyll concentrations. Similar results, albeit with slightly stronger responses of red SIF to variation in chlorophyll content, have been reported using the 1D SCOPE model (Rossini et al., 2016). To determine if these differences are due to parameter selection or due to differences in the underlying canopy structure (1D vs 3D) requires further investigation. In addition to chlorophyll, changes in foliar carotenoid content, dry matter content, and senescent factor also had an effect on the TOC SIF ( Fig. 11). Carotenoid contents had a proportionally similar negative effect on both red and far-red fluorescence, which can be explained in terms of competition for APAR with chlorophyll within the FLUSPECT model. Note however that certain carotenoids may also operate as accessory pigments passing excitation energy to chlorophyll, which was not considered here. For dry matter and senescent factor, the negative impact was stronger for far-red SIF compared to red SIF, which, in addition to the competition for APAR at red wavelengths, would point out towards differences in the specific absorption coefficients used in FLUSPECT for these two factors between 685 and 740 nm (Vilfan et al., 2016). In contrast, leaf water content did not directly affect TOC SIF (Fig. 11c) due to its lack of absorption in the fluorescence emission region (Jacquemoud et al., 2009). Indirect effects of leaf water content on SIF, which could be expected via physiological controlling mechanisms (Flexas and Medrano, 2002;Wohlfahrt et al., 2018), were not considered here.
Scaling up, the sensitivity analysis highlighted a strong influence of canopy structure parameters on the TOC SIF signal. All fluorescence emission wavelengths were sensitive to changes in overstory and understory LAI, although the magnitude of the response of SIF to changes in LAI was slightly greater at far-red relative to red wavelengths (Fig. 10c). This result is also in agreement with previous studies based on SCOPE modelling scenarios (Du et al., 2017). The contrasting response between wavelengths can be explained in terms of enhanced canopy re-absorption of red fluorescence at higher LAI values. In turn, the response of SIF to LAD was the same order of magnitude as that of both LAI and fCover (Fig. 10b). This previously published impact (Du et al., 2017;Migliavacca et al., 2017;Zeng et al., 2017) illustrates the importance of considering possible spatial variations in species-specific LAD (Pisek et al., 2013)   , suggested that red SIF is more sensitive to photosynthesis than far-red SIF. Although both red and far-red SIF were influenced by the overall quantum yield of fluorescence (or Fqe) (Fig. 11f), a parameter shown to respond to seasonal variation in photosynthesis (Ensminger et al., 2004;Soukupová et al., 2008;Springer et al., 2017), far-red fluorescence was more responsive to variations in LAI, FCover, understory LAI and especially foliar chlorophyll content (Fig. 13). These results suggest that far-red TOC SIF may better capture the full space of factors affecting photosynthesis (related to both APAR and Fqe), whereas red TOC SIF is perhaps less sensitive to variations in APAR, albeit with a greater sensitivity to Fqe. Considered together, and acknowledging the caveat that our results are based on a local sensitivity analysis conducted at a single forest stand, our results illustrate the wavelength-dependent action of multiple and scale-dependent structural, biochemical and physiological processes coexisting within a heterogeneous forest canopy and underlying the interpretability of TOC SIF.

Implications, limitations and next steps
We presented a data-driven tool, based on the DART model, capable of simulating reflectance (BRF) and SIF in an accurate virtual forest scene that accounted for a geometrically realistic description of branch, tree and forest scale architecture. Although being beyond the scope of this study, the next logical step is to validate the modelling results with high spatial resolution measurements of canopy SIF across different forest types. In this context, our simulated level of TOC far-red SIF (0.86 mW/m 2 /sr/nm), using the default empirical data for parameterization, was within the range of reported measurements in deciduous forest sites and at peak growing season: 0.8-1.0 mW/m 2 /sr/nm for a mixed forest of red oak (Quercus rubra) and yellow birch (Acer rubrum) (Yang et al., 2015), or average values of 0.4, 1.2, and 1.7 (mW/ m 2 /sr/nm) for stands dominated by linden, maple and oak, respectively (Rossini et al., 2016).
Despite its promise, the TLS-based reconstruction technique applied here has a few notable limitations. Firstly, owing to the resolution limit of terrestrial laser scanner, some tiny branches and twigs are too small to be detected by TLS and reconstructed. Additionally, self-occlusion hides some woody elements. Finally, as the leaf intersection algorithm used here is based on the reconstructed woody elements and an assumption that the leaves are uniformly distributed around the branches (Åkerblom et al., 2018), it may lead to an inaccurate leaf distribution within crown.
With regards to the limitations of the study, our sensitivity analysis did not consider multi-factorial influences. Due to computational constraints, it was limited to a (single variable at a time) local analysis. Such a limitation could result in a potentially inadequate understanding of the simulated TOC SIF signal. In a near future, a more robust global sensitivity analysis (e.g., Verrelst et al., 2015Verrelst et al., , 2016 should be feasible, given the projected increases in DART computational efficiency. Contrary to 1D SCOPE ( Van der Tol et al., 2009), the DART model does not include a leaf photosynthetic and energy balance modelling, due to its complexity in 3D space. Instead, we varied the Fqe parameter which, together with APAR, acts as the main control between SIF and photosynthesis. Possible coupling of DART with a leaf-level photosynthesis model would allow to take into account the spatial and temporal variation of leaf-level biochemical and physiological parameters within a crown, their response to local light environment and their impact in leaf and canopy-level photosynthetic rates. These applications could help to improve our understanding of how the TOC SIF signal is related to diurnal and seasonal changes in GPP, but also open up new possibilities to study relevant ecological questions such as the functional role of sunflecks (Pearcy, 1990;Way and Pearcy, 2012), or the coupling between canopy structural and functional traits (Chambers et al., 2007;Hardiman et al., 2013).

Conclusions
We presented a study modelling TOC SIF with consideration of the fine 3D structural heterogeneity found in a forest ecosystem. We do this by constructing 3D forest scene based on TLS data and inserting it in the DART canopy radiative transfer model, which has been coupled with the FLUSPECT leaf fluorescence model. Both models were parametrized using empirical data measured for the forest stand of interest. The potential of this scheme was demonstrated with a local sensitivity analysis, which revealed new insight into the factors that control TOC SIF in a complex forest ecosystem. We conclude that understory SIF can contribute significantly to TOC SIF in open ecosystems, which should be considered when relating SIF and GPP. We also found evidence that, contrary to expectations, red TOC SIF may be a more sensitive indicator of leaf-level Fqe variations due to its smaller dependency on the actual chlorophyll content, a factor which determines the compensation point between amount of APAR and SIF reabsorption effects.