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 Univercity, Sapporo, Japan 3Global Station for Food, Land, and Water Resources, GI-Core, Hokkaido University, Sapporo, Japan Correspondence: Yuma Sakai (ysakai@jamstec.go.jp)


Introduction
Global terrestrial ecosystems control the atmospheric CO 2 concentration through gross primary production (GPP) and ecosystem respiration processes (Canadell et al., 2007;Richardson et al., 2009;Piao et al., 2013). The ecosystem responses to climate 25 change have not yet been adequately quantified because of insufficient observations and modeling ability (Bunn and Goetz, 2006;Lasslop et al., 2010). Thus, there is great demand in the scientific community for methods of constraining global GPP through existing observation networks (Anav et al., 2015;Teubner et al., 2019). Estimating GPP is essential for various applications, ranging from yield predictions to evaluating and predicting the impact of regional and global environmental changes (Waring et al., 1998;Schimel, 2007). 30 Chlorophyll fluorescence is an energy release pathway for excess incident light in the photosynthetic process. Over the last ten years, extensive studies have revealed that canopy-scale sun-induced chlorophyll fluorescence (SIF) can be observed from satellites, such as the Greenhouse gases Observation Satellite (GOSAT-1&2) (Frankenberg et al., 2011), Orbiting Carbon Observatory-2 (OCO-2) (Li et al., 2018;Norton et al., 2019), Global Ozone Monitoring Experiment-2 (GOME-2) (Joiner et al., 2013), and TROPOspheric Monitoring Instrument (TROPOMI) (Köhler et al., 2018) using Fraunhofer lines in the near-35 infrared spectral domain. Satellite-derived SIF potentially provides a direct pathway linking leaf-level photosynthesis to global GPP (Guanter et al., 2010;Frankenberg et al., 2011;Joiner et al., 2013;Porcar-Castell et al., 2014). For example, the observed SIF exhibits a good correlation with net photosynthesis, which is quantified by the monitoring gas exchange method at the leaf level and using the eddy covariance method at the ecosystem scale (Wieneke et al., 2018). SIF can be used to infer the photosynthetic capacity of the plant canopy (Zhang et al., 2018). However, it is not clear how leaf-level SIF emissions 40 contribute to the top-of-canopy directional SIF, because satellite-observed SIF uses the near-infrared spectral domain, in which multiple scattering on the leaf surface is dominant. Based on the steady-state fluorescence yield theory (Genty et al., 1989), a model for leaf-level SIF and photosynthesis under various environmental conditions has been developed ( Van der Tol et al., 2014). The spectral variability of emitted SIF radiance has also been quantified by a radiative transfer model at the leaf level (Pedrós et al., 2010), canopy level (Tol et al., 2009;Gastellu-Etchegorry et al., 2017;Yang and van der Tol, 2018;Liu et al., without sufficient input data, it is difficult to extend the model simulations to the various ecosystems around the world. This paper describes a 3-D 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.
2 Model description 70 We developed a 3-D 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)). As one of the series of FLiES modules, FLiES-SIF is a radiative transfer model for solar radiation in the visible and near-infrared domains (Kobayashi and Iwabuchi, 2008) and for thermal emissions in the thermal infrared domain (Kobayashi et al., 2012). The FLiES-SIF model shares several of the key aspects of numerical schemes in FLiES: it employs a spatially explicit forest landscape and is 75 based on a Monte Carlo ray-tracing approach. Multiple scatterings among leaves, trunks, and soil background are numerically simulated with an unbiased approach (Kobayashi and Iwabuchi, 2008). The simulated landscapes are represented by spatially explicit geometric tree crown objects (see details in Sect. 2.4). The performance and reliability of FLiES for simulating light transmittance through a canopy and bidirectional reflectance factors have been 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 SIF radiance at the top of the 80 canopy, but the simulation codes can easily be extended to simulate SIF at any height level within the plant canopy.

Bidirectional SIF radiance
The bidirectional SIF radiance at wavelength λ at the top of the atmosphere, I (λ, Ω V ), can be decomposed into four different light transfer pathways: I (λ, Ω V ) = I dir_sun (λ, Ω V ) + I dir_shade (λ, Ω V ) + I ms_sun (λ, Ω V ) + I ms_shade (λ, Ω V ) (1) 85 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. (1) indicate direct SIF radiance from sunlit leaves, direct SIF radiance from 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 90 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 mentioned, the wavelength λ is omitted from subsequent equations.
The intensity of SIF is related to the absorbed photosynthetically active radiation (APAR c ) taken in by the forest canopy.

95
If the forest is sparse or the leaf area density in the tree crowns is low, a large portion of incident PAR is attenuated through 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 arbitrary APAR 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. The actual APAR can be independently calculated for a given canopy landscape. In subsequent sections, we describe the radiance components derived with 105 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 4π 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. 120

Attenuation function
The attenuation of SIF in the view direction Ω V is calculated by the attenuation function exp (−τ V ). When the hotspot effect is negligibly small, the attenuation function is expressed using the plant canopy gap fraction theory: where u, s, G, and γ are the leaf area density, pathlength, mean leaf projection area, and clumping index, respectively. The 125 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 scale. Because FLiES-SIF version 1.0 employs explicit tree crown landscapes, clumping larger than the crown scale need not 130 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 in a previous report (Kobayashi et al., 2010). The hotspot effect refers to the strong illumination near the solar direction

135
(Ω 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
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 Tol et al., 2014). In FLiES-SIF 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 150 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) 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 160 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.4.

Phase function for SIF emissions
The phase function for SIF emissions P f gives the fraction of SIF emitted in the view direction Ω V . Similar to the scattering 165 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.

Multiple scattering
SIF emissions from the leaf surface occur in all directions (upward and downward in the plant canopy), although they are not always isotropic, as shown in Sect. 2.2.3. A certain portion of SIF does not directly go toward the sky. This portion hits other 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 will return to the view direction, which enhances the observed SIF radiance depending on the magnitude of the multiple-scattering 175 contribution. The multiple-scattering process of SIF is the same as the scattering process of solar radiation, and the multiplescattering component can be formulated in exactly the same way as the bidirectional reflectance factor described in Kobayashi and Iwabuchi (2008). The 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 (14) has exactly the same form as Eq. (16) in Kobayashi and Iwabuchi (2008). The form of the phase function P (Ω j , Ω V ) is also described by Eq. (7) in Kobayashi and Iwabuchi (2008).

185
The attenuation function is the same as described in Sect. 2.2.1.

Canopy structure represented by FLiES-SIF version 1.0
The simulated landscapes are represented by the spatially explicit geometric tree crown objects (Fig. 1). Three crown objects (spheroid, cone, and cylinder) are defined for the SIF simulations. These crown objects are further separated into leafy-crown and woody domains: the outer and inner domains are filled with leaves and wood, respectively. In the default setting, the height 190 and diameter of the woody domain is set to be half that of the outer domain ( Fig. 1). Stems are represented by solid cylinders.
The individual tree dimensions can be defined differently. The default landscape size is 100 m × 100 m, though it is possible to change these values. To create the virtual forest canopy for SIF simulations, it is necessary to determine all of the tree positions in the forest. If there are ground-based tree census data, they can be used to create the virtual forest canopy. The virtual forest canopies are reconstructed by a statistical approach . Assuming that the spatial distribution of trees follows a

195
Poisson or Neyman distribution, the individual tree positions are determined by these statistical functions and random numbers.
FLiES has a module for the voxel representation of the forest landscape Wu et al. (2018). However, this module is currently not incorporated into FLiES-SIF version 1.0.

Photon tracing algorithm
The numerical scheme of the photon tracing is shown in Fig. 2 In FLiES-SIF, 3-D forest landscapes are reconstructed using geometric tree objects composed of cones, cylinders, and spheroids (see Sect. 2.5 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 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 210 leafy-canopy voxels 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 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 1m 3 . Note that the extracted leafy voxels are not always completely filled with canopy geometry because of the presence of canopy edge voxels, which only partially contain the leafy-215 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.

B. 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.
The position v 0 is determined as follows. First, an arbitrary voxel is chosen at random from the voxel table (Fig. 3). The exact 220 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 225 contains branch domains, the randomly determined position v 0 may be outside the leafy canopy. Therefore, the position v 0 is 8 https://doi.org/10.5194/gmd-2020-19 Preprint. Discussion started: 21 February 2020 c Author(s) 2020. CC BY 4.0 License. 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.

C. Determination of the leaf properties for SIF emission
After position v 0 has been determined, the leaf properties at the selected position are determined. Two leaf properties are 230 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. 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. (20), the program returns to Eq. (19) and calculates another θ L . In Eq. (20), 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.

D. Compute the leaf-level SIF emission and the direct SIF radiance in the view direction
Once the position v 0 and leaf properties have been determined, the leaf-level SIF emission w 0 and the direct SIF radiance 250 (I dir_sun and I dir_shade ) can be computed using the equations derived in Sects. 2.1 and 2.2. For the sunlit leaf condition, the calculation of w 0 includes the spherical integration of the attenuation function (Eqs. (10) and (11)), which is time-consuming.

E. Determination of the new emission direction
Direct SIF radiance in the view direction Ω V is determined by procedure D. 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 265
Test simulations of the SIF emissions were performed on a one-hectare virtual forest (Fig. 4). The aim of these tests was to understand the sensitivity of the SIF simulated by FLiES-SIF with respect to the geometric conditions (solar zenith angle, SZA; view zenith angle, VZA), sunlit leaf fraction, and LAI, and to identify the factors (hotspots, light attenuation, phase function, weight of photons) that contribute to SIF radiance under the given forest structure. The individual tree positions and sizes were determined at random. The spheroid shape was employed for the individual crowns. The tree density used 270 in the sensitivity analysis was 359 trees ha −1 . The canopy layer height was set to 25 m (Fig. 4) 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. The model requires optical data in the PAR domain and the spectral wavelength to be simulated. In this sensitivity analysis, we used the data assembled by Kobayashi (2015a). Figure 5 shows the spectral leaf reflectance and transmittance and the woody/soil reflectance. The leaf reflectance and transmittance, woody reflectance, and soil reflectance were calculated 275 from various broadleaf spectral data, medium reflective woody elements, and medium reflective soil surfaces in Kobayashi (2015b), respectively. All optical data were averaged over 10-nm intervals between 650 nm and 850 nm. The optical data in the PAR domain were computed as the average from 400-700 nm ( Table 1) (Fig. 7) and Farquhar et al. (1980). Tol's model is based on energy partition within leaves. In calculating φ f , the photosynthetic capability is measured by applying the pulse amplitude modulated fluorometry (PAM) system to actual leaves.

290
However, we obtained the photosynthetic rate in the sensitivity analysis using Farquhar's model (Farquhar et al., 1980) instead of the PAM measurement process, because our analysis considers a virtual forest. The parameter values in these models were set by reference to previous literature (such as Van der Tol et al. (2014) and De Pury and Farquhar (1997)), and the results compared with those using a constant APAR-dependent (Tol's model) φ f . The incident total PAR on the canopy surface was fixed at 2000 µ mol m −2 s −1 , except for in the APAR sensitivity analysis, and the fraction of diffuse radiation was fixed at 0.3.

295
In the sensitivity analysis, we used 10 5 -10 6 photons in each model run. We conducted three transition tests to study the model response to changes in LAI (0-20), VZA (−80°-80°), and SZA (0°-75°). Figures 8 and 9 indicate the dependency of SZA and LAI on total SIF radiance in each wavelength. In the following section, we analyze the sensitivity of SZA and LAI in more detail using the results for λ=760 nm.
3.1 Angular dependency of SIF 300 Figure 8 shows that the total SIF radiance for wavelengths between 650-840 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 (e.g., Fig. 8(a)), 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 direction of forward emission (VZA > 0), SIF increases with an increase in VZA and sharp strong peaks 305 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.
To analyze the dependency of SZA in more detail, we explored the influence of SZA on three terms in Eq. (1), namely the 310 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 10 shows the dependence on SZA 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 domains are minor. In contrast, I ms_sun + I ms_shade exhibits bowl-like shapes (Fig. 10(d)), which contribute to the enhancement 315 of total SIF at higher angles. In the FLiES-SIF model framework, SIF radiance is computed by collecting the contribution factor (Eqs. (4) and (14)) 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. (1)-(9)) 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. (3.3.1)).
When α j is sufficiently large and the hotspot effect is marginal, the attenuation function is determined by the forest structure 320 (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 ) contribute to the total SIF increase in higher angular domains. In addition, I dir_shade makes a slightly larger 325 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. 11. 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.1), and then adjusted by multiplying by the ratio of APAR c to APAR app (Eq. (2)).
The simulated angular patterns indicate that APAR c increases with an increase in SZA ( Fig. 12(b)). The increase in APAR c 335 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. 12(a)). This is because APAR c 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 photon path (L p ) increases, the gap fraction P sun becomes smaller (Eq. (16)). As a result, shaded leaves are more likely to 340 be selected in the random process in Eq. (17). 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. (1), which means the statistical accuracy of I dir_sun decreases as SZA increases. Depending on the target sampling variables to be simulated, the number of photons should be determined (i.e., more 345 photons may be necessary to investigate the behavior of I dir_sun in cases where the sunlit leaf fraction is low). Figure 9 shows the sensitivity of total SIF radiance to LAI for wavelengths of 650-840 nm. The simulated SIF increases with LAI and then becomes saturated over the whole wavelength range, although the speed of saturation varies with the wavelength. In the visible domain ( Fig. 9(a)-(d)), the simulated SIF becomes saturated when LAI = 2. In the near-infrared domain ( Fig. 9(k)-(y)), 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. (1), namely I dir_sun , I dir_shade , and I ms_sun + I ms_shade , as well as the total SIF radiance (I) (forward direction in Fig. 13 and backward direction in Fig. 14). In our simulation scenarios, I dir_sun contributed about 54% of total SIF radiance when LAI = 3, VZA = 10°, and SZA = 20°. I dir_hade and I ms_sun + I ms_shade contributed 7% and 39%, respectively (Fig. 15). Figures 13 and 14 show 355 that the individual SIF components respond differently to the LAI.

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. 13(b) and 14(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 360 (Fig. 16(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. ) has a major influence on simulated direct SIF radiance from sunlit leaves. The SIF radiance increases and then becomes saturated without decreasing when α j is equal to 0, because the rate of decrease in I dir_sun becomes small when τ = 0. Additionally, 365 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 of the increase in the pathlength to the canopy boundary when passing through other canopies. The influences of |VZA| and LAI are prominent 370 in negative VZA directions. In this case, the decrease in SIF radiance with an increase in LAI becomes significant because of the SIF emitted through the local canopy to the view point, and the attenuation in the local canopy (and inn other canopies) increases with LAI. Thus, the increase in pathlength as |VZA| increases significantly affects I dir_sun in the view direction.

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 de-375 creasing when VZA is negative (Figs. 13(c) and 14(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 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 380 range where ψ rapidly decreases with any increase in τ .

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 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. 13(d) and 14(d)). When VZA is positive, the SIF radiance increases and then becomes 385 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.2), 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 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 390 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 expectation of the product of w and exp(−τ ) is larger than when VZA is negative (Fig. 14). 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 395 a similar reason as for I dir_shade when SZA is negative. Figure 16 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 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 400 visible SIF radiance ( Fig. 9(a)-(d)). 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.
3.4 Influence of fluorescence yield on variable APAR L scenario Figure 17 compares the total SIF derived from fixed and variable leaf-level SIF yields. To explore the influence of the SIF 405 yield on the above dependencies, we derive φ f by means of Tol's model (Van der Tol et al., 2014) to calculate the yield from APAR, which is obtained by our model (Fig. 7). Additionally, we used Farquhar's model (Farquhar et al., 1980) to obtain data on the photosynthesis rate using PAM observations ( Van der Tol et al., 2014). Figure 17 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.

410
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.
In this paper, behind the SIF emissions from complex forest canopies. We performed a test run to demonstrate the sensitivity of SIF to the view angle, LAI, and leaf-level SIF yield. The simulation results show that SIF increases with LAI before becoming saturated when LAI > 2-4, depending on the spectral wavelength. The sensitivity analyses also showed that simulated SIF radiation may decrease with LAI when LAI > 5. These phenomena were observed under certain sun and view angle conditions. This type of nonlinear and nonmonotonic SIF behavior with respect to LAI is also related to the spatial forest structure pattern.

425
The hotspot effect plays an important role in SIF simulations when the view direction is close to the sun direction. The SIF yield φ f influences the canopy SIF, especially when APAR is low. In FLiES-SIF version 1.0, the leaf-level SIF yield model is not directly coupled: the SIF yield should be determined from the literature or existing models for use as an input variable.
FLiES-SIF version 1.0 can be used to quantify the canopy SIF at various view angles, including the contribution of multiple scatterings, which is an important component in the near-infrared domain. The proposed model can be used to standardize 430 satellite SIF by correcting the bidirectional effect. This step will contribute to improved GPP estimation accuracy through SIF.
In this model description paper, we have focused on the formulation and simulation schemes of FLiES-SIF version 1.0, and have presented the results from sensitivity analyses of major variables such as LAI. Model validation using field measurements will be performed in future studies. Thorough validation against measured quantities should be conducted to evaluate the accuracy of the model.    (d) Figure 13. 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 θS = 20°and φS = 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 14. 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 θS = 20°and φS = 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.