FLiES-SIF ver. 1.0: Three-dimensional radiative transfer model for estimating solar induced fluorescence

estimating solar induced fluorescence Yuma Sakai1,2, Hideki Kobayashi1, and Tomomichi Kato2,3 1Institute of Arctic Climate and Environment Research, Research Institute for Global Change, Japan Agency for Marine-Earth Science and Technology, Yokohama, Japan 2Research Faculty of Agriculture, Hokkaido University, Sapporo, Japan 3Global Station for Food, Land, and Water Resources, GI-Core, Hokkaido University, Sapporo, Japan Correspondence: Yuma Sakai (ysakai@jamstec.go.jp)

and Energy fluxes (SCOPE) model, Tol et al. (2009)) have been widely used to analyze the physiological, meteorological, and geometrical influences on observed SIF. These plane parallel models provide some insight into the general mechanisms 60 behind the temporal and seasonal variations in SIF. However, the lack of complexity in their actual canopy structures means that 1D models often give inaccurate directional SIF features. Three-dimensional (3D) models (Zarco-Tejada et al., 2013;Gastellu-Etchegorry et al., 2017;Hernández-Clemente et al., 2017), although requiring vast computational resources, have the potential to delineate the realistic directional canopy SIF. The radiative transfer model used in SIF simulations should exhibit several characteristics. First, the contribution of sunlit and shaded leaves to canopy-scale directional SIF emissions should be separately quantified. The intensity of SIF depends on the absorbed photosynthetically active radiation (APAR) on leaf surfaces, and the emissions from sunlit and shaded leaves are quite different (the APAR of sunlit leaves can be 100 times higher than that of shaded leaves). Second, the multiple scattering of fluorescence should be accurately computed, as most satellites use the near-infrared spectral domain. Third, although 3D models are required to evaluate realistic SIF features, the model's input variables should be easily created or accessible from existing databases. This is because, without sufficient input data, 70 it is difficult to extend the model simulations to the various ecosystems around the world. This paper describes a 3D Monte Carlo plant canopy radiative transfer model, the Forest Light Environmental Simulator (FLiES) for simulating canopy-scale directional SIF radiance, and evaluates the performance of the model by analyzing the angular and multiple-scattering effects on SIF. We developed a 3D plant canopy radiative transfer model for simulating the canopy-scale directional SIF radiance (Forest Light Environmental Simulator for SIF, FLiES-SIF version 1.0, Kobayashi and Sakai (2019)). Figure 1(a) shows the overall framework of the FLiES family modules. The aim of FLiES is to consider the impact of landscape heterogeneity on the 80 radiative processes and determine how this links to the canopy energy, water, and carbon exchanges. One of the important aspects of modeling is to make the modules as computationally efficient as possible: most of the FLiES modules, including the newly proposed FLiES-SIF, can be easily run on moderate personal computers with relatively modest memory demands and public software (GNU gfortran, gcc, and R). The model development was initiated using 3D radiative transfer modeling in the shortwave domains (Kobayashi and Iwabuchi, 2008). A 1D version of the atmospheric radiative transfer model, MCARaTS 85 (Iwabuchi, 2006), was incorporated to simulate the atmosphere-forest radiation interaction. A longwave radiative transfer module was then added, together with energy balance and plant physiology modules ( Fig. 1(a), Kobayashi et al. (2012); Baldocchi and Harley (1995)). All these modules are related to the radiation emitted in the Stefan-Boltzmann law of the sun, by the earth's surface, and by atmospheric media.
The current FLiES-SIF work adds a radiation interaction module for the induced radiation emitted from leaf pigments, and describes how to combine the forest structure information and leaf physiology models (Van der Farquhar et al., 1980) with the fluorescence radiative transfer module ( Fig. 1(a)). The FLiES-SIF model shares some key aspects of numerical schemes in FLiES: it employs a spatially explicit forest landscape (Sect. 2.1.2) and is based on a Monte Carlo raytracing approach (Sects. 2.2-2.3). Analogous to the modeling in previous FLiES modules, FLiES-SIF employs a Monte Carlo sampling scheme, where photon-tracing sequences represent the integration form of the radiative transfer equation, or the so-95 called Neumann's series (Antyufeev and Marshak, 1990). In such modeling, the photon path lengths and scattered directions are determined by random numbers and probability distribution functions such as the Lambert-Beer exponential function and the scattering phase function. The scattering and re-absorption of emitted fluorescence light must also be considered to identify the relationship between the fluorescence emitted by the chloroplasts and the top-of-canopy outgoing fluorescence radiance (Porcar-Castell et al., 2014). Several recent studies have worked on the quantification of the impact and modeling of 100 scattering and absorption effects from the leaf scale (e.g., Agati et al. (1993); van der Tol et al. (2019)) to the canopy scale (e.g., Romero et al. (2018)). Multiple scatterings and re-absorption among leaves, trunks, and soil background can be numerically simulated using unbiased and efficient approaches (Kobayashi and Iwabuchi, 2008). The performance and reliability of FLiES for simulating light transmittance through a canopy and bidirectional reflectance factors have been extensively investigated in previous studies (Widlowski et al., 2011(Widlowski et al., , 2013(Widlowski et al., , 2015. As a default setting, FLiES-SIF version 1.0 simulates the bidirectional 105 SIF radiance at the top of the canopy, but the simulation codes can easily be extended to simulate SIF at any height level within the plant canopy and for the spatial mapping of SIF radiance at the top of the canopy.

Canopy structure represented by FLiES-SIF
The forest landscapes employed in FLiES-SIF consist of simple geometric objects such as cones, spheroids, and cylinders, as in other 3D models (e.g., DART and FluorFLIGHT) (Fig. 2). The volume domains inside the crown objects can be further split 110 into multiple domains to realize the spatial distribution of the leaf area and woody area densities. This approach is useful in some aspects because (1) it can establish realistic plant canopies in which the majority of leaves are distributed in the outer and upper parts of the crowns (see Fig. 2), and (2) this approach is simple and computationally efficient. The insides of the crown volumes are assumed to be turbid media, where the light attenuation follows the Lambert-Beer exponential law as a function of leaf/woody area densities and leaf inclination angles. The conventional turbid medium approach assumes that the 115 leaves are randomly distributed in space. In the FLiES modules, including FLiES-SIF, a spatially anisotropic arrangement of leaves (the so-called clumping effect) is modeled using re-collision probability theory Stenberg, 2003, 2005;Kobayashi et al., 2010), which is particularly important for the shoot-scale clumping of needle leaf. More details can be found in the FLiES user manual (Kobayashi, 2019). Note that FLiES has a module for the voxel representation of the forest landscape (Wu et al., 2018), but this module is not currently incorporated into FLiES-SIF version 1.0.

120
The realization of individual tree objects has some degree of freedom in the characterization of leaf and woody densities in crowns. In the sensitivity analysis described in the next chapter, we consider the following forest architectures. The crown objects are separated into two domains, namely the outer leafy-crown and inner woody domains; the outer and inner domains are filled with 100 % leaves and 100 % wood, respectively. The height and diameter of the inner woody domain is set to be half that of the outer domain (Fig. 2). Stems are represented by solid cylinders. The individual tree dimensions can be defined 125 differently. The landscape size used in this study is 100 m × 100 m. To create the virtual forest canopy for SIF simulations, it is necessary to determine all of the tree positions in the forest. If ground-based tree census data exist, they can be used to create the virtual forest canopy. The virtual forest canopies can also be reconstructed using a statistical approach .
Assuming that the spatial arrangement of the trees follows a Poisson or Neyman distribution, the individual tree positions are determined by these statistical functions and random numbers. The simulation flow of the spectral SIF calculation is shown in Fig. 1(b). The FLiES-SIF model requires four major inputs, namely geometry data, meteorological data, forest stand data, and optical data for leaves and other elements (e.g., soil background). The geometry data specify the sun and sensor view directions (zenith and azimuth angles). The meteorological data are used in the precomputation of fluorescence yield. The incident total and diffuse photosynthetically active radiation (PAR) 135 data from the top of forest canopies are also used in the canopy radiative transfer module. If no PAR observations are available, these data can be calculated by the FLiES atmosphere module (1D MCARaTS, see Fig. 1(a)).
Before running the Monte Carlo ray-tracing processes in the forest canopy, the forest structures (leaf area density and leaf voxel look-up table) and leaf-level physiology (fluorescence module) are computed. In FLiES-SIF, landscape-scale LAI is an input variable. The model requires the leaf area density of individual canopy volumes. FLiES-SIF computes the leaf area 140 density from a given forest landscape and LAI data. Fluorescence ray-tracing also requires information about the leaf-level sun-induced fluorescence yield and its spectral composition. This information is computed by the leaf photosynthesis and fluorescence module. These modules are described in the following subsections (Sects. 2.1.4 and 2.1.5).
Once the forest structure and leaf fluorescence yield have been computed, Monte Carlo ray tracing is performed in the broad PAR domain to determine the PAR absorbed in the forest landscape (APAR c in Sect. 2.2). The SIF radiance is then simulated 145 on an individual wavelength basis. Details of the radiative transfer formulation and ray-tracing algorithm are summarized in the following sections (Sects. 2.2 -2.5).

Creation of the leafy-canopy voxel look-up table
In FLiES-SIF, 3D forest landscapes are reconstructed using the geometric tree objects described in Sect. 2.1.2 and Fig. 3 upper). The photon tracing starts from an arbitrary position v 0 = (x, y, z) within the leafy-canopy volume. This position is 150 determined by three random numbers corresponding to x, y, and z. When the canopy landscape is sparse, the majority of randomly determined positions v 0 will be outside of the leafy-canopy space, which means a large number of trial runs will be required to determine an appropriate position v 0 . To reduce the computation time, regularly placed leafy-canopy voxel tables are extracted to determine where to start the SIF emission and subsequent photon tracing (Fig. 3). In FLiES-SIF version 1.0, the leafy canopy voxel information is saved in a look-up table (Fig. 3). The voxel information in the table contains lower and upper   155 corner positions (x, y, z and x + dx, y + dy, z + dz), the leaf area density (LAD) of the voxel, and the sunlit leaf area density (LAD sun ). The size of each voxel is currently 1m 3 . Note that the extracted leafy voxels are not always completely filled with canopy geometry-canopy edge voxels only partially contain the leafy-canopy geometry. In addition, tree canopy geometries contain branch domains. Thus, even if the voxel is completely inside the canopy geometry, there may be some domains that do not contain leaves.

Computation of leaf-level fluorescence yield
In FLiES-SIF version 1.0, the leaf-level sun-induced fluorescence yield (hereafter SIF yield φ f ) is pre-computed and stored in a look-up table prior to ray tracing. The SIF yield φ f is computed by the model of Van der Tol et al. (2014) and Farquhar's photosynthesis model (Farquhar et al., 1980) under various environmental conditions, including absorbed PAR (APAR). The actual photosynthesis can not be determined without stomata models (e.g., Collatz et al. (1991)). The leaf temperature is also 165 dependent on photosynthesis and stomata regulations. These interrelations are solved by iterating the energy balance, stomata, and photosynthesis equations. The CANOAK-FLiES module ( Fig. 1(a), Kobayashi et al. (2012)) can handle such leaf-level coupled physiology phenomena, but this would require more input variables and increase the computational load. Thus, in the current FLiES-SIF, we adapted the following assumptions to obtain reasonable photosynthesis simulation results. First, the leaf temperature was assumed to be the same as the surface air temperature. This is usually acceptable, except in very 170 dry conditions when the stomata are almost closed in daytime. The other assumption concerns the stomata modeling. The FLiES-SIF module does not explicitly use the stomata model. Rather, the consequences of the stomata activity, i.e., the downregulation of the intercellular partial CO 2 pressure (ipCO 2 ), were modeled as a function of the vapor pressure deficit (VPD).
We used the experimental relationships measured by Dang et al. (1997), who investigated the relationships between ipCO 2 and VPD in three tree species (pine, spruce, and aspen, see Fig. 10 of Dang et al. (1997)). If we simulate SIF over such species, 175 the regression lines given by Dang et al. (1997) can be used. For other species, we created a simplified function based on the relationship derived by Dang et al. (1997): where apCO 2 denotes the ambient partial CO 2 pressure.
To simulate the spectral SIF, the spectral composition of SIF must be known. Our approach is similar to that used in the SCOPE model (Tol et al., 2009). We derived the spectral composition from the FluorMODleaf model (Zarco-Tejada et al., 2006;Pedrós et al., 2010). The calculated leaf-level spectral SIF radiance variations given by FluorMODleaf were normalized to determine the fraction of SIF at wavelength λ, f s mWm −2 sr −1 , with respect to the broadband SIF Wm −2 . The 185 standard-setting and variables are described in Sect. 3.1. That is, we only used the fraction of spectral composition from the FluorMODleaf model. The radiance was then determined from APAR and φ f , which varies with environmental conditions and leaf traits such as the maximum carboxylation capacity, V cmax , used in the photosynthesis model.

Bidirectional SIF radiance
The bidirectional SIF radiance at wavelength λ at the top of canopy, I (λ, Ω V ), can be decomposed into four different light 190 transfer pathways: where the subscripts "dir" and "ms" indicate the direct emission of SIF and SIF after multiple scatterings, respectively. The direction vector Ω V = (θ V , φ V ) contains the observation zenith and azimuth angles. The radiance elements I dir_sun , I dir_shade , I ms_sun , and I ms_shade on the right-hand side of Eq.
(2) indicate direct SIF radiance from sunlit leaves, direct SIF radiance from 195 shaded leaves, sunlit SIF radiance after multiple scatterings, and shaded SIF radiance after multiple scatterings, respectively.
Here, the direct emission of SIF indicates SIF that is emitted from leaves and directly escapes from the canopy space without hitting other leaves and trunks. On the contrary, "multiple scattering SIF" indicates SIF that is emitted from leaves, hits other leaves, trunks, or soil background, and then escapes from the canopy space in the view direction. Note that most of the optical and radiance quantities described below are spectral variables. For simplicity of the mathematical expressions, if not explicitly 200 mentioned, the wavelength λ is omitted from subsequent equations.
The intensity of SIF is related to the absorbed photosynthetically active radiation (APAR) taken in by the forest canopy. If the forest is sparse or the leaf area density in the tree crowns is low, a large portion of incident PAR is transmitted through plant canopies. The transmitted PAR does not contribute to the SIF emissions on the leaf surface. Thus, if photon tracing is performed under sparsely vegetated canopies, the simulation includes large amounts of photons that are not used to compute SIF. To make 205 the numerical simulation more efficient, we propose a variance reduction technique. FLiES-SIF forces all incident PAR to be absorbed by sunlit or shaded leaves and initiates the photon tracing for SIF emitted from leaves. This procedure artificially enhances or diminishes APAR, biasing the simulated SIF depending on the ratio of actual APAR to the "apparent APAR" (APAR app ) used in the simulation. Thus, the simulated SIF under the APAR app is adjusted to the actual APAR (APAR c ) conditions: where I (Ω V ) and I (Ω V ) denote the SIF radiance with APAR c and APAR app , respectively. APAR app is simulated with the SIF simultaneously. The APAR c is independently calculated for a given canopy landscape before the SIF simulation ( Fig. 1(b) and Sect. 2.1.3). In subsequent sections, we describe the radiance components derived with APAR app (I dir_sun , I dir_shade , I ms_sun , and I ms_shade ).

Calculation of direct SIF radiance
The direct SIF radiance from sunlit and shaded leaves is calculated by summing all direct SIF radiation contribution factors of the i-th photon (ψ dir,i ): where V sun , V shade , v 0 , and N indicate the classes of sunlit and shaded leaves, the position of the photon (x, y, z), and the total number of photons, respectively.
The direct SIF radiation contribution factor of the i-th photon ψ dir, i can be decomposed into three components: leaf-level SIF emission weight w 0 , directional emission transfer function (the so-called phase function P f ), and attenuation function: Here, τ V is the optical thickness of the plant canopy in the view direction Ω V . The factor 2π is a normalization factor for the phase function P f . These three components in ψ dir , namely w 0 , P f , and exp (−τ V ), indicate the SIF emitted in all directions from both adaxial and abaxial sides of a single leaf, the fraction of SIF emitted in the view direction, and the fraction of SIF attenuation to the top of the canopy in the view direction, respectively.

Attenuation function 230
The attenuation of SIF in the view direction Ω σ is calculated by the attenuation function exp (−τ V ). When the hotspot effect is not considered, the attenuation function is expressed using the plant canopy gap fraction theory: where u i , s i , G σ,i and γ i are the leaf area density, path length, mean leaf projection area, and clumping index of the i-th tree. They are aggregated over the trees located in the light path between the emission point to the top of canopy in the view 235 direction, respectively. The path length, s, is a sum of canopy paths that penetrates through crown objects.
The mean leaf projection area G is a function of the leaf inclination angle distribution function g L and an arbitrary direction Ω σ (such as the sun direction Ω S or view direction Ω V ): Generally, the clumping index contains various nonrandom scales of spatial leaf distributions, from the shoot to the landscape 240 scale. Because FLiES-SIF version 1.0 employs explicit tree crown landscapes, clumping larger than the crown scale need not be considered. However, the crown volumes are expressed as turbid media: if the leaves are not randomly distributed in the crown object, e.g., shoot-scale clumping (Cescatti and Zorer, 2003;Chen et al., 1997), attenuation must be corrected according to the shoot-scale clumping index. In FLiES-SIF version 1.0, the shoot-scale clumping index is estimated by the spherically averaged shoot silhouette area (Cescatti and Zorer, 2003). Details on how shoot-scale clumping is incorporated can be found 245 in a previous report (Kobayashi et al., 2010). The hotspot effect refers to the strong illumination near the solar direction (Ω V ≈ −Ω S ). When the hotspot effect is nonnegligible, the modified optical thickness τ is expressed as: where H is a hotspot function expressed by the Hapke model (Hapke, 2012), which is used in the framework of the FLiES model (Kobayashi and Iwabuchi, 2008): where Ω j , l, and α j indicate the incident direction after the j-th scattering, the radius of the disk-shaped flat leaves, and the scattering angle (α j = cos −1 |Ω V · Ω j |), respectively.

Leaf-level SIF emission weight 255
The leaf-level SIF emission weight w 0 can be calculated from the SIF yield φ f and APAR on the leaf surface (APAR L ): where f s is the fraction of SIF at wavelength λ (mW m −2 sr −1 ) with respect to the broadband SIF (W m −2 ). Thus, f s is a function of wavelength. The SIF yield φ f is a function of APAR L and various environmental and leaf trait variables such as ambient air temperature, humidity, CO 2 concentration, and carboxylation capacity (Van der . In FLiES-SIF 260 version 1.0, φ f is read from a look-up table across a wide range of APAR L , which should be pre-computed by the leaf-level SIF yield models. The exact computation of APAR L under the angular dependency of PAR can be performed by backward ray tracing at the given position of a leaf, but this approach is time-consuming. For more efficient simulations, the values of APAR L for sunlit and shaded leaves are approximated as the product of the incident-diffuse PAR and the attenuation function exp (−τ s) 265 integrated over the upper hemisphere: where PAR dir and PAR dif denote the incident direct and diffuse PAR, respectively, ω PAR is the average single-scattering albedo in the PAR spectral domain (400 -700 nm), and ω PAR is the sum of the leaf reflectance r PAR and transmittance t PAR in the PAR domain (ω PAR = r PAR + t PAR ). This equation assumes that diffuse PAR is isotropic over the sky and neglects direct PAR scattered within the plant canopy and soil background. Thus, APAR L may be underestimated when the background reflectance is high, such as in the case of snow cover. To further reduce the computation time, the hemispherical integration of the attenuation function is approximated by an average of the limited-angle samplings. Details of the computation method are given in Sect. 2.5.

Phase function for SIF emissions 275
The phase function for SIF emissions P f gives the fraction of SIF emitted in the view direction Ω V . Similar to the scattering phase function for the reflection of solar illumination, P f can be determined by the following equations: where f ada and f aba are the fraction of SIF emissions from the adaxial and abaxial sides of a leaf; f ada + f aba = 1. Note that, in our definition, we have assumed that illumination by solar beams is always on the adaxial side of a leaf. leaves, trunks, or soil background. The SIF energy from those impacts is scattered, goes in another direction, and then impacts something else. We define this process as the multiple scatterings of SIF. After multiple scatterings, some of the SIF energy 285 will return to the view direction, which enhances the observed SIF radiance depending on the magnitude of the multiplescattering contribution. The multiple-scattering process of SIF is the same as the scattering process of solar radiation, and the multiple-scattering component can be formulated in exactly the same way as the bidirectional reflectance factor described in Kobayashi and Iwabuchi (2008). The scattered SIF radiance emitted by sunlit and shaded leaves is defined as I ms_sun and I ms_shade , respectively, and these radiance contributions can be calculated by summing all of the scattering contributions: Here, ψ i,j is calculated as follows: where w i,j is the weight of the i-th photon after the j-th scattering obtained by using the single-scattering albedo in the SIF spectral domain ω SIF = r SIF +t SIF (w i,j = w i,j−1 ω SIF ). Equation (15) is exactly the same as the multiple scatterins in the shortwave 295 radiative transfer (Kobayashi and Iwabuchi, 2008). The form of the phase function P (Ω j , Ω V ) is also described by Eq. (7) in Kobayashi and Iwabuchi (2008). The attenuation function is the same as described in Sect. 2.3.1.

Photon tracing algorithm
The numerical scheme of the photon tracing is shown in Fig. 4. The procedures framed by the dotted grey rectangle indicate the photon tracing scheme for direct SIF emissions. The area outside the dotted grey rectangle corresponds to scattered photon 300 tracing. The algorithm for scattered photon tracing is exactly the same as the photon tracing method for solar radiation. Here, we focus on the SIF emission scheme in the grey rectangle. Details of the scattered components are summarized in Kobayashi and Iwabuchi (2008).

A. Set a new photon in the leafy-canopy
The position v 0 = (x, y, z) from which SIF emission occurs within a leafy-canopy domain is determined by random numbers.

305
The position v 0 is determined as follows. First, an arbitrary voxel is chosen at random from the voxel table (Fig. 3). The exact position (x, y, z) within a selected voxel is then determined by three random numbers (Rx, Ry and Rz; R ∈ [0, 1]): where v l = (x l , yl, zl) denotes the position of the lower corner of the selected voxel. If the selected voxel is an edge voxel or contains branch domains, the randomly determined position v 0 may be outside the leafy canopy. Therefore, the position v 0 is checked to determine whether it is in the leafy domain. If the position is outside the leafy domain, the program generates a new random number and selects another voxel. This procedure continues until the leafy canopy position v 0 is obtained.

B. Determination of the leaf properties for SIF emission 315
After position v 0 has been determined, the leaf properties at the selected position are determined. Two leaf properties are required to continue the computation of the SIF emission: the leaf illumination status (sunlit or shaded) and the leaf surface normal vector Ω L = (θ L , φ L ). The sunlit leaf area fraction P sun at v 0 is computed using the interception of direct sunlight: where L p is the cumulative LAI at v 0 along the path of the sunlight and G S is a mean leaf projection area defined in Eq. 7. The leaf illumination status (sunlit or shaded) is then determined by a random number R: The leaf surface normal vector Ω L is also required because the leaf-level SIF emission is related to APAR at the leaf surface (APAR L ). APAR L is computed from the cosine of the sunlight and leaf normal angles. Assuming the leaves are randomly distributed, the azimuthal angle of the leaf surface normal φ L can be determined by: For a given leaf angle distribution function g L := g (θ L ), the zenith angle of the leaf surface normal θ L can be determined by the rejection method. In the first step, θ L is calculated using a random number: Then, θ L is further evaluated using g L : If θ L is rejected by the abovementioned criteria in Eq. (21), the program returns to Eq. (20) and calculates another θ L . In Eq. (21), the evaluation function is a form of leaf angle distribution function multiplied by a sine value. This sine comes from the Jacobian of the polar coordinate and is necessary because g L is defined in polar coordinates.
C. Compute the leaf-level SIF emission and the direct SIF radiance in the view direction 335 Once the position v 0 and leaf properties have been determined, the leaf-level SIF emission w 0 and the direct SIF radiance (I dir_sun and I dir_shade ) can be computed using the equations derived in Sects. 2.2 and 2.3. The calculation of w0 includes the spherical integration of the attenuation function (Eqs. (11) and (12)), which is time-consuming. Thus, FLiES-SIF version 1.0 approximates this spherical integration by taking the average of five directions (θ, φ) = (0°, 0°), (60°, 0°), (60°, 90°), (60°, 180°), and (60°, 270°). We tested the performance of this 5-angle assumption by comparing with 10-degree interval samplings 340 (9 zenith and 36 azimuth angles = 324 angle sampilngs). When the attenuation functions were computed by these two angle samplings at 104 randomly selected positions in the forest landscapes used in the sensitivity analysis in section 3, the mean absolute error of this approximation was 14.6 % (N = 10000). Finally, I dir_sun and I dir_shade are calculated by the local estimation method using Eqs. (4) and (5) (Antyufeev and Marshak, 1990;Marchuk et al., 1980).

D. Determination of the new emission direction 345
Direct SIF radiance in the view direction Ω V is determined by procedure C. The multiple scattering contribution is further evaluated by photon tracing. To start the photon tracing, the emission direction Ω (θ, φ) is calculated using two random numbers and the leaf surface normal vector Ω L = (θ L , φ L ). Assuming that the SIF emission is bi-Lambertian on the leaf surface, the zenith and azimuthal angles relative to the leaf normal (α, β) are determined by: The scattering direction Ω (θ, φ) in the Cartesian coordinate system is then calculated by a coordinate transformation from (α, β) to (θ, φ).

Sensitivity analysis
We created a one-hectare virtual forest as a default conditions (Fig. 5 (a)). Fourthly, we performed 1D and 3D comparisons in an actual diurnal variations in total and diffuse PAR observed in Yokohama, Japan (35°22'N, 139°37'E) in the summer of 2014 (Delta-T sunshine sensor, Delta-T Co Ltd.) (Sect. 3.5). In this exercise, we evaluated the potential errors (overestimation) in the 1D homogeneous layer (turbid medium) approach. We compared 365 in seven different scenarios that include different leaf angle distributions (spherical, erectrophile, and planophie) and withincrown clumping (γ = 0.6) in the 3D landscapes (Table 1). Lastly, we evaluated the sensitivity of the leaf level fluorescence yield (Sect. 3.6). In this exercise, we computed the leaf level fluorescence with the FLiES-SIF physiology module ( Fig. 1(b)).

Input data and simulation condition
The individual tree positions and sizes were determined at random. The spheroid shape was employed for the individual 370 crowns. The tree density used in the sensitivity analysis was 359 trees ha −1 . The canopy layer height was set to 25 m (Fig. 5) and the crown coverage was 96 %. FLiES-SIF assumes that all crowns have the same leaf area density. The spherical leaf angle distribution function was used. We also used three different landscape conditions, tropical broad leaf ( Fig. forest (b)), Evergreen needleleaf ( Fig. 5 (c)), and Savanna ( Fig. 5 (d)) to understand the effect of landscape structure on SIF. The model requires optical data in the PAR domain and the spectral wavelength to be simulated. In this sensitivity analysis, we used the data 375 assembled by Kobayashi (2015a). Figure 6 shows the spectral leaf reflectance and transmittance and the woody/soil reflectance.
The leaf reflectance and transmittance, woody reflectance, and soil reflectance were calculated from various broadleaf spectral data, medium reflective woody elements, and medium reflective soil surfaces in Kobayashi (2015b), respectively. All optical data were averaged over 10nm intervals between 650 nm and 850 nm. The optical data in the PAR domain were computed as the average from 400 -700 nm ( Table 2). The same woody reflectance data were used for both stem surface and branch 380 materials. The fractions of SIF emission (f ada , f aba ) were determined using the FluorMODleaf model (FluorMODgui V3.1) (Zarco-Tejada et al., 2006;Pedrós et al., 2010) (Fig. 7). To run FluorMODgui V3.1, we used the default biochemical parameters (leaf structure parameter N = 1.5, chlorophyll a+b content C ab = 33.0 µ g cm −1 , water content C w = 0.025 cm, dry matter content C m = 0.01 g cm −2 , fluorescence quantum efficiency F i = 0.01, leaf temperature T = 20.0 • C, species temperature dependence = 2 (beans), stoichiometry of PSII (photosystem II) to PSI reaction centers Sto = 2.0) under the downward spectral 385 sky radiation data (direct transmittance in sun direction (τ S ), FluorMOD30V23.MEP). The fractions of SIF emission were derived from the simulated leaf fluorescence output by normalizing the simulated leaf level SIF from the adaxial and abaxial sides. In this sensitivity analysis, we employed two types of leaf-level SIF yield φ f . The first type is a constant value of 0.01 throughout the whole APAR range from Sect. 3.2 -3.5. This value is used to test the impact of forest structures (LAI) and sun and observation geometries on SIF. The second type is an APAR-dependent value derived using the models of Van der Tol et al.

390
(2014) (Fig. 8) and Farquhar et al. (1980 analysis, and the fraction of diffuse radiation was fixed at 0.3. In the sensitivity analysis, we used 10 5 -10 6 photons in each model run. Figure 9 indicates the dependency of SZA and LAI on total SIF radiance between 650 to 850 nm. In the following section, we analyze the SIF sensitivity in λ = 760 nm.

Intercomparisons with the DART model
To demonstrate the efficacy of FLiES-SIF for calculating the SIF radiation, we compared the output with that of DART 400 (Gastellu-Etchegorry et al., 2017), which is one of the available 3D models. We adopted the flux-tracking mode in the radiative method of DART ver. 5.7.6 to calculate the SIF radiation, and compared the simulation results under the following conditions: SZA = 30°, SAA (Solar Azimuth Angle) = 0°, and λ = 760 nm. Only part of the default landscape (40 m × 40 m, 50 trees) was simulated to reduce the computational load. Figure 10 compares the SIF radiation calculated by FLiES-SIF with that given by DART-SIF. In terms of VZA dependency, FLiES-SIF overestimates SIF radiation by about 18 % on average in 405 the forward direction (VZA > 0), and the difference becomes larger as VZA increases. In contrast, FLiES-SIF underestimates the SIF radiation by about 12 % in the backward direction (VZA < 0), and the difference reaches a maximum at VZA = −50°.
In terms of LAI dependency, FLiES-SIF overestimates SIF radiation by about 9 % on average, with the difference being especially pronounced when LAI > 6. The SIF radiation increases with LAI in FLiES-SIF, but decreases as LAI increases in the DART model. As a result, FLiES-SIF gives similar SIF radiation values as DART, and has greater sensitivity to angles.

410
FLiES-SIF is a reasonable and proper model for determining the SIF radiation. The DART model has a useful GUI and can calculate SIF radiation on complex landscape structures using many kinds of 3D objects. However, unlike FLiES-SIF, DART does not include a leaf physiological module. Additionally, FLiES-SIF requires less computational resources than the DART model.

Sensitivity analysis with default forest landscape 415
3.3.1 Angular dependency of SIF Figure 9(a) shows the total SIF radiance for wavelengths between 650 -850 nm. These figures indicate that the SIF radiance shows a strong peak near the sun direction over the whole wavelength range, although the SZA value, which exhibits the maximum SIF, varies according to the wavelength. In the visible red region, the SIF radiance reaches an extremum at a lower SZA than in the near-infrared region. Regardless of wavelength, the angular dependency of SIF exhibits similar patterns: in the 420 direction of forward emission (VZA > 0), SIF increases with an increase in VZA and sharp strong peaks appear around the sun direction (the hotspot effect). In the backward direction (VZA < 0), the SIF decreases with an increase in |VZA| and attains a minimum at around −35°-−50°, before increasing with |VZA|. Although the general angular patterns are similar across the whole wavelength range, the strength of the hotspot peak in the forward direction and the minimum SIF in the backward direction vary slightly with the wavelength.

425
To analyze the dependency of VZA in more detail, we explored the influence of VZA on three terms in Eq.
(2), namely the direct SIF from the radiance of sunlit and shaded leaves (I dir_sun , I dir_shade ) and the scattered radiance (I ms_sun + I ms_shade ), as well as the total SIF radiance (I). The simulated SIF shows distinct angular features for each SIF component (I dir_sun , I dir_shade , I ms_sun + I ms_shade ). Figure 11 shows the dependence on VZA of the SIF components when LAI = 3.0 for a wavelength of 760 nm. I dir_sun has a strong peak near the sun direction because of the hotspot effect, whereas angular changes of I dir_sun in other 430 domains are minor. In contrast, I ms_sun + I ms_shade exhibits bowl-like shapes (Fig. 11(d)), which contributes to the enhancement of total SIF at higher angles. In the FLiES-SIF model framework, SIF radiance is computed by collecting the contribution factor (Eqs. (5) and (15)) from the attenuation function, weight of photons, and phase function. Among those factors, the drastic changes in the optical thickness of the attenuation function (Eqs. (2) -(10)) contributed the most to the hotspot in I sun_dir . The attenuation function displays a strong peak around the sun direction because of the hotspot parameter (H in Eq. (9)). When α j 435 is sufficiently large and the hotspot effect is marginal, the attenuation function is determined by the forest structure (such as LAI and leaf angle density). Away from the sun direction, the SIF radiance gradually decreases or increases slightly. This angular feature (VZA) is influenced by the initial photon weight and phase function through the dependency on the leaf surface normal: the initial photon weight is calculated as the inner product between the leaf angle and the sun direction. The influence of SZA on the phase function is greater than that on the initial photon weight. The other two components (I dir_shade , I ms_sun + I ms_shade ) 440 contribute to the total SIF increase in higher angular domains. In addition, I dir_shade makes a slightly larger contribution in the backward direction, because shaded leaves tend to be more aligned with the backward direction. The shaded leaves only absorb diffuse sky radiation, so the relative magnitude of I dir_shade with respect to I dir_sun greatly depends on the fraction of diffuse radiation. The contributions of these three components to the direct in four difference sun angles are presented in Fig. 12. These partitions vary with the fraction of incoming diffuse radiation, optical properties (leaf reflectance and transmittance, woody and soil reflectance), and the leaf area.

Angular dependencies of APAR and sunlit leaves
Because SIF radiance is greatly affected by the APAR of the leaves, the angular behavior of APAR is essential in understanding the numerical computation of SIF emissions. In the FLiES-SIF model, the SIF radiance is first computed under the apparent APAR (APAR app ) conditions (Sect. 2.2), and then adjusted by multiplying by the ratio of APAR c to APAR app (Eq. (3)).

450
The simulated angular patterns indicate that APAR c increases with an increase in SZA ( Fig. 13(b)). The increase in APAR c with respect to SZA corresponds to the increase in the photon pathlength inside the forest canopy. As SZA increases, more photons are likely to hit leaves before they pass through the canopy layers. In contrast, APAR app decreases as SZA decreases ( Fig. 13(a)). This is because APAR app is related to the fraction of sunlit leaves. As described in simulation procedure C in Sect. 2.5, the photon tracing is initiated from either sunlit or shaded leaves at randomly selected positions. As the LAI along the 455 photon path (L p ) increases, the gap fraction P sun becomes smaller (Eq. (17)). As a result, shaded leaves are more likely to be selected in the random process in Eq. (18). In other words, as the fraction of shaded leaves increases, the amount of energy in the simulated system decreases. In the Monte Carlo simulations, the statistical accuracy of the simulated variables depends on the number of photon samplings. The decrease in APAR app does not affect the simulated accuracy of the total SIF radiance; however, it does affect the individual components in Eq.
(2), which means the statistical accuracy of I dir_sun decreases as SZA 460 increases. Depending on the target sampling variables to be simulated, the number of photons should be determined (i.e., more photons may be necessary to investigate the behavior of I dir_sun in cases where the sunlit leaf fraction is low).

Leaf area density dependency
Figure 9(b) shows the sensitivity of total SIF radiance to LAI for wavelengths of 650 -850 nm. The simulated SIF increases with LAI and then becomes saturated over the whole wavelength range, although the speed of saturation varies with the 465 wavelength. In the visible domain, the simulated SIF becomes saturated when LAI = 2. In the near-infrared domain, the simulated SIF is not saturated at higher LAI values, indicating that SIF is more sensitive to LAI in the near-infrared domain.
To analyze the dependency on LAI, we explored the influence of LAI on three terms in Eq.

470
I dir_hade and I ms_sun + I ms_shade contributed 7 % and 39 %, respectively (Fig. 16). Figures 14 and 15 show that the individual SIF components respond differently to the LAI.

i. Direct SIF radiance from sunlit leaves
The LAI dependency of direct SIF radiance from sunlit leaves is influenced by the hotspot function and the magnitude of VZA (Figs. 14(b) and 15(b)). Generally, the SIF radiance emitted from sunlit leaves increases and then saturates as LAI increases, because the number of sunlit leaves also increases and becomes saturated, although the fraction of sunlit leaves decreases ( Fig. 17(c)). However, in terms of simulated SIF radiance, there are ranges of LAI in which SIF radiance decreases with an increase in LAI. In these regions, the decrease in SIF radiance is caused by the attenuation of SIF radiance in the canopy.
The magnitude of this attenuation depends on both the hotspot function and VZA. The hotspot function (i.e., the angle α j in Eq. 9) has a major influence on simulated direct SIF radiance from sunlit leaves. The SIF radiance increases and then 480 becomes saturated without decreasing when α j is equal to 0, because the rate of decrease in I dir_sun becomes small when τ = 0.
Additionally, smaller values of α j produce a smaller rate of decrease in I dir_sun with respect to increases in LAI through the hotspot effect. The magnitude of VZA (i.e., |VZA|) also influences the simulated SIF radiance. Generally, larger LAI values lead to a decrease in the attenuation of SIF radiation from sunlit leaves in the canopy when VZA is positive, because most sunlit leaves inhabit the canopy surface. However, the attenuation of SIF radiation in other canopies increases with |VZA| because  ii. Direct SIF radiance from shaded leaves The fraction of shaded leaves has a major influence on SIF radiance. SIF increases and then becomes saturated without decreasing when VZA is negative (Figs. 14(c) and 15(c)). This variation in SIF is caused by an increase in the fraction of shaded leaves, because the rate of increase in the fraction is larger than the rate of decrease in I dir_shade . In contrast, the rate of decrease in I dir_shade becomes greater than the rate of increase in the fraction of shaded leaves when VZA is positive. In this region, the 495 expectation of the pathlength to the view point is larger than for negative VZA, because the canopy surface is covered with sunlit leaves. This increase in optical thickness, which depends on the pathlength, has a major effect on I dir_shade in the LAI range where ψ rapidly decreases with any increase in τ .
iii. Scattered SIF radiance The scattered SIF radiance refers to the sum of the scattered radiance from sunlit and shaded leaves, I ms = I ms_sun + I ms_shade 500 in our model. The LAI dependency with respect to view direction on the scattered SIF radiance is in contrast to the direct radiance from shaded leaves (Figs. 14(d) and 15(d)). When VZA is positive, the SIF radiance increases and then becomes saturated without decreasing. The pathlength from sunlit leaves to the population boundary in the view direction has a major influence on simulated scattered SIF radiance. As previously explained (Sect. 3.3.3), the surface of the canopy is covered by sunlit leaves, which provide a large photon weight to scattered photons, in the positive VZA direction. When LAI is large, the 505 decrease in I ms with an increase in LAI becomes vanishingly small. This is because the scattered radiation from high-weight photons reaches the view point with little attenuation. Larger values of LAI lead to shorter scattering pathlengths and fewer scatterings, so the photon weight w j is larger. Additionally, the pathlength between sunlit leaves and the boundary of the canopy is nearly constant, irrespective of LAI variation.
The simulated SIF radiance, therefore, becomes larger than the radiance in the negative VZA direction. Actually, the expec-510 tation of the product of w and exp(−τ ) is larger than when VZA is negative (Fig. 15). In contrast, with an increase in LAI, the SIF radiance decreases and becomes saturated after increasing because of the increase in τ from sunlit leaves. This is for a similar reason as for I dir_shade when SZA is negative. Figure 17 shows APAR and the fraction of sunlit leaves as a function of LAI. APAR c increases with an increase in LAI and becomes saturated at around LAI = 2. APAR app and the fraction of sunlit leaves decrease when LAI < 2. The increase 515 in APAR c and the decrease in APAR app are more abrupt than the SIF increase with respect to LAI. This is because APAR is the visible light where the absorption of green leaves is high (˜0.9). Thus, the APAR c saturation curve has similar patterns of visible SIF radiance. At higher LAI, the fraction of sunlit leaves is low and APAR app decreases. The statistical accuracy of I dir_sun becomes drastically lower as APAR app and the fraction of sunlit leaves decrease. Accurate simulations of I dir_sun require an increased number of photons to be traced. 520

Variation of landscape
FLiES-SIF is applicable to a wide variety of landscapes. Figure 18 shows the LAI and SZA dependencies of the total SIF radiation in three different landscapes, namely tropical broadleaf forest, evergreen needle forest, and savanna, which have different tree densities, shapes, and crown coverages . These simulations use same optical data to demonstrate the applicability of our model. The SIF radiation increases with LAI and then becomes saturated, and there is little difference 525 in LAI dependency among these landscapes (upper panels in Fig. 18). In contrast, the angular dependency is affected by the landscape composition. For example, in the case of the savanna (Fig. 18(c)), the hotspot peak is smaller than for the other landscape types because the attenuation within other crowns is reduced by the lower crown coverage.

Comparison of 1D and 3D actual diurnal variations in PAR
Simulations with a one-hour step size were performed using observed total and diffuse PAR data from Yokohama, Japan from 530 July 12 -18 , 2014 ( Fig. 19(a)). This period included clear sky days (DOY 192,195,and 196) and overcast days (DOY 193,198). The LAI value was fixed to 3.0 throughout the simulation period in all scenarios. The diurnal variations in sunlit LAI show distinct features with respect to 1D and 3D landscapes, as well as to the leaf angle distribution (Fig. 19(b)).
In summary, the simulated sunlit LAIs in the 1D scenarios (1DsphNC, 1DereNC, 1DplaNC) are higher than those in the 3D scenarios (3DsphNC, 3DereNC, 3DplaNC, 3DsphWC), and the difference in the leaf angle distribution affects the diurnal 535 feature of sunlit LAI significantly: erectrophile cases give the highest values, followed by spherical and then planophile. The planophile cases (1DplaNC, 3DplaNC) exhibit weak diurnal variations because of the high probability of horizontal leaves: the uppermost leaves receive most of the incoming light, regardless of SZA. The sunlit LAI of the within-crown clumping scenario (3DsphWC) is significantly lower than in the no-clumping scenario (3DsphNC, 29.7 % lower at noon). Figures 19(c) and 19(d) show the hourly simulation results for the top-of-canopy SIF at 760 nm.

540
The simulated SIF generally follows the diurnal variations in incoming PAR; however, the absolute range varies greatly depending on the leaf angle distribution (spherical, erectrophile, or planophile) and in the 1D and 3D cases. For example, the canopy SIF simulated by 1DplaNC is 2.3 times higher than that of 3DereNC at noon on a clear day (DOY 192). In contrast, the simulated SIF differences are minor on overcast days (DOY 193,198). The simulated results also indicate that there is a sort of "trade-off" relationship between the amount of sunlit LAI and SIF: the SIF values simulated by planophile scenarios 545 (1DplaNC, 3DplaNC) are higher than those in erectrophile scenarios (1DereNC, 3DereNC), while the amount of sunlit LAI in planophile scenarios is lower than that in erectrophile scenarios. This is because leaves with an erectrophile distribution receive the solar beam efficiently as they have a large proportion of vertical leaves. This, in turn, results in larger incident angles of the solar beam on the leaf surfaces, which reduces the incident PAR intensity on the unit leaf area. Comparing the 3D scenarios with the 1D scenarios, the latter simulated the SIF to be 33 -41 % higher than the former.

550
Additionally, we investigated the contribution of direct SIF from sunlit leaves (I dir_sun ), shaded leaves (I dir_shade ), and multiple scattering components (I ms_sun + I ms_shade ) (Figs. 19(e) and 19(f)). The largest contribution comes from sunlit leaves. The contribution from shaded leaves is generally lower than that from sunlit leaves, although on overcast days (DOY 193 and 198) the contribution becomes close to or even higher than that of the sunlit leaves (Figs. 19(e) and 19(f)). Multiple scattering contributes substantially in the near-infrared domain (30 -40 % of total SIF radiance), but is expected to make a lower contribution 555 in the red spectrum because of low leaf reflectance and transmittance. A comparison of 3DsphNC and 3DsphWC shows that the shaded leaves of the within-crown clumping scenario contribute more to the total SIF than in the no-clumping scenario. Overall, the 1D and 3D comparisons, as well as differences in leaf angle distributions, suggest that reasonable assumptions must be made regarding the canopy structure if the SIF simulations are to be reliable. If 1D models are applied to forest canopies, the simulated SIF is prone to be overestimated. 3.6 Influence of fluorescence yield on variable APAR L scenario Figure 20 compares the total SIF derived from fixed and variable leaf-level SIF yields. To explore the influence of the SIF yield on the above dependencies, we derive φ f by means of Tol's model (Van der  to calculate the yield from APAR, which is obtained by our model (Fig. 8). Additionally, we obtain the photosynthesis yield by Farquhar's model using parameter values set by reference of previous literatures (such as Van der Tol et al. (2014) and De Pury and Farquhar (1997)) 565 to derive the φ f without observation data. Figure 20 compares the total SIF radiance I between models based on constant φ f and APAR-dependent φ f in terms of their dependence on LAI and SZA, respectively. The dependency on the two parameters is not substantially different because the variation in φ f is smaller than that of APAR. However, φ f affects APAR app as well as the bidirectional SIF radiance. Thus, obtaining accurate values of φ f is important in estimating the exact level of SIF; this issue will be considered in future work.   Springer, 1980. Norton, A. J., Rayner, P. J., Koffi, E. N., Scholze, M., Silver, J. D., and Wang, Y.-P.: Estimating global gross primary productivity using geometry data, meteorological data, forest stand data, and optical and leaf trait data, are required to run the model.         20°, (c) 40°, and (d) 60°. The contribution of shaded leaves is basically small and the contribution rates of the other two radiances exhibit some angular dependency. In the backward direction, the contribution of scattered radiation to SIF is greater than in the forward direction. (d) Figure 14. LAI dependency on SIF radiance. (a) Total radiance I, (b) direct radiance from sunlit leaves Idir_sun, (c) direct radiance from shaded leaves Idir_shade, and (d) radiance after multiple scatterings Ims_sun + Ims_shade at SZA = 20°and SAA = 0°(forward direction). Each line represents a different VZA value (0°-70°). SIF radiance increases with LAI and then becomes saturated in most cases. However, when VZA is large (e.g., black and red lines), the direct radiance decreases with an increase in LAI. (d) Figure 15. LAI dependence of SIF radiance. (a) Total radiance I, (b) direct radiance from sunlit leaves, (c) direct radiance from shaded leaves, and (d) radiance after multiple scatterings at SZA = 20°and SAA = 180°(backward direction). Each line represents a different VZA value (0°-70°). SIF radiance increases with LAI and then becomes saturated in most cases. However, when VZA is large (e.g., black and red lines), only the direct radiance from sunlit leaves decreases with an increase in LAI, different from the forward direction case.  (right)). The contribution of shaded leaves is small and the contribution of scattered radiance increases with LAI and VZA. Additionally, in the backward direction, the contribution of scattered radiation to SIF is larger than in the forward direction, similar to the angular dependency.