Simulation-Based Evaluation of the Estimation Methods of Far-Red Solar-Induced Chlorophyll Fluorescence Escape Probability in Discontinuous Forest Canopies

The escape probability of Solar-induced chlorophyll fluorescence (SIF) can be remotely estimated using reflectance measurements based on spectral invariants theory. This can then be used to correct the effects of canopy structure on canopy-leaving SIF. However, the feasibility of these estimation methods is untested in heterogeneous vegetation such as the discontinuous forest canopy layer under evaluation here. In this study, the Discrete Anisotropic Radiative Transfer (DART) model is used to simulate canopy-leaving SIF, canopy total emitted SIF, canopy interceptance, and the fraction of absorbed photosynthetically active radiation (fAPAR) in order to evaluate the estimation methods of SIF escape probability in discontinuous forest canopies. Our simulation results show that the normalized difference vegetation index (NDVI) can be used to partly eliminate the effects of background reflectance on the estimation of SIF escape probability in most cases, but fails to produce accurate estimations if the background is partly or totally covered by vegetation. We also found that SIF escape probabilities estimated at a high solar zenith angle have better estimation accuracy than those estimated at a lower solar zenith angle. Our results show that additional errors will be introduced to the estimation of SIF escape probability with the use of satellite products, especially when the product of leaf area index (LAI) and clumping index (CI) was underestimated. In other results, fAPAR has comparable estimation accuracy of SIF escape probability when compared to canopy interceptance. Additionally, fAPAR for the entire canopy has better estimation accuracy of SIF escape probability than fPAR for leaf only in sparse forest canopies. These results help us to better understand the current estimation results of SIF escape probability based on spectral invariants theory, and to improve its estimation accuracy in discontinuous forest canopies.


Introduction
Solar-induced chlorophyll fluorescence (SIF) is emitted from within the photosynthetic apparatus of all higher plants and is a good indicator of carbon assimilation and plant physiological status [1,2]. in structure and contain overstory and understory layers with varying species and phenologies which depend on factors such as site fertility [27].
In this study, we investigate whether the estimation methods of far-red SIF escape probability based on SIT holds for discontinuous forest canopies using the model-based assessment strategy. To achieve this, we used a 3-D RT model, DART, to simulate the TOC SIF and per-leaf emitted SIF (i.e., total emitted SIF) based on the 3-D radiative budget simulation results. After that, the escape probability of SIF in the forest area was directly calculated based on the 3-D simulation results. Subsequently, the performance of estimation methods for escape probability based on SIT was analyzed under different vegetation conditions, such as different leaf optical properties, canopy structure, and background optical properties.

Spectral Invariants Theory and SIF Escape Probability
The original estimation method (Equation (1)) for SIF escape probability at far-red band based on SIT was developed by [22]. It combines top of canopy reflectance NIR v , canopy interceptance i o , and leaf single scattering albedo ω λ as In a remote sensing application, NIR v is generally replaced by NIR T or NIR T ·NDVI, where NDVI is used to solve the "black soil problem" [12]. In addition, fAPAR can be used to replace canopy interceptance i o because fAPAR is simpler to estimate. Based on these considerations, we decided to calculate and evaluate SIF escape probability ( f esc SIT : f esc based on SIT) using four formulas shown in Table 1.

Simulation of SIF Escape Probability
The Discrete Anisotropic Radiative Transfer (DART) model is a 3-D RT model [28] that works with realistic representations of plant architecture described by the foliage distribution, including leaf clumping at branch and crown levels, along with detailed geometry of stems and branches. Recently, the leaf fluorescence FLUSPECT model has been integrated into the DART model to simulate fluorescence from leaf to canopy level [18]. Besides the simulation of spectral total and SIF radiance that escape the canopy, DART also simulates the 3-D, 2-D, and 1-D spectral radiative budgets of vegetation scenes per type of scene element (e.g., ground, woody elements, and leaves). It can also simulate the radiative budget of individual facets (i.e., triangles) that are used to represent scene elements such as leaves. These radiative budgets are computed at different stages of the radiative transfer modeling: after direct sun illumination, after atmosphere illumination, and after each subsequent DART iteration. They include radiation that is intercepted, scattered, absorbed, and emitted (i.e., SIF in the short waves and thermal emission in the long waves), and also the irradiance and exitance at each face of each voxel with which the landscape is simulated. Units can be W/m 2 /µm or percentage of top of canopy irradiance. Therefore, DART can compute the light regime and instantaneous incoming photosynthetically active radiation (PAR) and fAPAR of each leaf element and each voxel and layer of the canopy.
The spectral SIF escape probability (f esc ) was defined as the ratio of canopy-leaving SIF observed by sensors (SIF obs ) to total emitted SIF from whole-canopy (SIF tot ) at a given wavelength [12], where SIF obs is the TOC SIF exitance. It corresponds to the angular integration of DART TOC SIF radiance images; more details about the simulation process can be found in [29]. SIF tot is the sum of all SIF photons emitted by all leaves within the canopy in all directions [12,24]. SIF tot was derived from the DART 3-D radiative budgets of SIF emitted in all directions for each leaf within the 3-D vegetation canopy.

Simulation of Canopy Interceptance
Canopy interceptance (i o ) is defined as the probability that an incoming solar photon interacts in first-order with leaves (or needles) within the canopy at any wavelength [21]. Because it is equal to one minus directional gap fraction, it is often derived from the simulation of directional gap fraction as in the works of [12,22] using the SCOPE model. In DART, i o is the "Intercepted" term of the 3-D and 1-D vegetation radiative budgets at the end of the direct sun illumination stage (i.e., "Illudir" stage in DART), before the scattering stage. It is equal to 1-"Intercepted" terms of the ground and trunk/branch radiative budgets of the Illudir stage.
where RB intercepted lea f and RB intercepted ground are the leaf absorbed irradiance and ground intercepted radiation at a specific wavelength, respectively.

Simulation of fAPAR
fAPAR canopy was defined as the fraction of absorbed photosynthetically active radiation (APAR) by the entire canopy (i.e., leaves and wood elements without ground) over the PAR region from 400 nm to 700 nm [30]. Accordingly, fAPAR leaf was defined as the fraction of absorbed photosynthetically active radiation (APAR) only by leaves. To evaluate the impact of woody elements on the simulation and estimation of SIF escape probabilities, two types of fAPAR (i.e., fAPAR canopy = fAPAR leaf+woody and fAPAR leaf ) were used to estimate SIF escape probabilities (Section 3.4).
In DART, fAPAR canopy and fAPAR leaf were derived from the radiative budget of facets per type of scene element (i.e., leaf, wood element, ground).
where E(λ) is the incident irradiance, and λ is wavelength over the spectral region from 400 to 700 nm, RB absorbed (λ) is the irradiance absorbed by leaf only (i.e., RB absorbed lea f (λ)) or entire canopy (i.e., RB absorbed canopy (λ) = RB absorbed lea f (λ) + RB absorbed woody (λ)). In addition, RB absorbed (λ) was calculated by area-weighted average as below: where n is the number of leaves in the canopy, RB absorbed i and S element i are the absorbed irradiance and area of scene elements of type i (e.g., leaf, woody element), respectively, and S element canopy is the sum of areas of all elements within the canopy.

Simulation Experiment
We conducted a numerical simulation experiment to simulate SIF escape probability using the 3-D DART model and then using it to evaluate the validity of the remote estimation methods described in Section 2.1. DART provides a realistic representation of plant architecture [18,29]. Hence, DART provides a way to simulate and evaluate SIF escape probability in a more realistic vegetation scene relative to simpler models, especially for sparse forest stands which are found across the boreal forest [31].
A 3-D Silver birch tree (Betula pendula Roth) reconstructed from Terrestrial laser scanning (TLS) data, measured in the vicinity of the Station for Measuring Ecosystem-Atmosphere Relations II (SMEAR II) in southern Finland (24.31478 • E, 61.84335 • N), was used to generate a 3-D simulation scene. We used an option of repetitive scene of the DART model to eliminate boundary effects [29]. The scene horizontal dimension was set to 4 m for default configuration (i.e., fractional vegetation cover (FVC) = 0.6, see Table 2). By adjusting the scene horizontal dimensions, the fractional vegetation cover (FVC) was modified from sparse forest (FVC = 0.3) to a dense forest (FVC = 0.8). Due to the computational constraints, four key input parameters of the DART model were varied to generate the simulation database (see Table 2), including leaf chlorophyll content (Cab), fractional vegetation cover (FVC), background spectra, and solar zenith angle (SZA). Variable values of leaf chlorophyll content were used to represent the variation of leaf optical properties. Varying values of fractional vegetation cover were used to describe the canopy structure of different forest types: 0.3, 0.6, and 0.8 are referred to as very sparse, little sparse, and dense forest, respectively. Given the two-layer structure of the forest and its seasonal variation in spectral properties, three typical background spectra were adapted to represent the change of background spectral properties. In addition, a non-reflecting background was used to test the ability of SIF escape probability that was estimated using SIT. Contrary to the DART simulations carried out by previous studies in which the solar zenith angle was set to a constant value [12,23], the solar zenith angle was varied here from 20 • to 70 • with a step of 10 • to represent the seasonal variation for most of the boreal forest areas. Details about the other default variables for the DART simulations can be found in [29].
All simulations were performed between 640 and 850 nm at a spectral resolution of 7.5 nm. Directional and total emitted SIF at 760 nm were used for the simulation and estimation of SIF escape probability. All canopy reflectance and SIF were retrieved in nadir viewing direction (i.e., view zenith angle equal to 0) to keep consistent with the default observation mode of current SIF satellite (such as the OCO-2 satellite). NDVI was derived from the canopy reflectance in the red (650 nm) and near-infra-red (NIR) (800 nm) bands according to [32]. SIF escape probabilities simulated by the DART Remote Sens. 2020, 12, 3962 6 of 17 model ( f esc DART ) were regarded as 'real' value, the relative estimation accuracy of SIF escape probabilities using SIT ( f esc SIT ) were calculated as ( f esc SIT − f esc DART )/ f esc DART .

Remote Estimation of Canopy Interceptance
In addition to canopy reflectance, canopy interceptance (i o ) is a crucial parameter for the estimation of SIF escape probability [22]. Although i o is not easily accessed directly from remote sensing measurements, it can be indirectly estimated from remotely sensed data using equations such that of [25]: where G(θ) is the mean value of leaf projection function for direction θ. It is generally fixed to 0.5 for global-scale remote sensing applications [33]. LAI and CI are the leaf area index and clumping index, respectively. They are directly available from satellite remote sensing products. To explore the uncertainties of the remote sensing products on the estimation of SIF escape probability, artificial errors (−50% and +50% constant value error of LAI*CI) were generated and then were used to estimate SIF escape probability.

Evalution Process
In Section 3.2, all four estimation formulas in Table 1 were evaluated with 'real' values simulated by the DART model (Section 3.1); in Section 3.3, only NIR v i• ·ω λ was used to evaluate the usability of remotely sensed canopy interceptance (i o ); in Section 3.4, only NIR v fAPAR·ω λ was used to evaluate the influence of woody elements on the estimation of SIF escape probability. The complete evaluation process is outlined as a flowchart in Figure 1.

Remote Estimation of Canopy Interceptance
In addition to canopy reflectance, canopy interceptance (io) is a crucial parameter for the estimation of SIF escape probability [22]. Although io is not easily accessed directly from remote sensing measurements, it can be indirectly estimated from remotely sensed data using equations such that of [25]: where ( ) is the mean value of leaf projection function for direction . It is generally fixed to 0.5 for global-scale remote sensing applications [33]. LAI and CI are the leaf area index and clumping index, respectively. They are directly available from satellite remote sensing products. To explore the uncertainties of the remote sensing products on the estimation of SIF escape probability, artificial errors (−50% and +50% constant value error of LAI*CI) were generated and then were used to estimate SIF escape probability.

Evalution Process
In Section 3.2, all four estimation formulas in Table 1 were evaluated with 'real' values simulated by the DART model (Section 3.1); in Section 3.3, only °• was used to evaluate the usability of remotely sensed canopy interceptance (io); in Section 3.4, only • was used to evaluate the influence of woody elements on the estimation of SIF escape probability. The complete evaluation process is outlined as a flowchart in Figure 1.

Simulation of 'Real' f esc in Discontinuous Forest Canopies
Using the variable input parameters listed in Table 2, TOC SIF and total emitted SIF were simulated respectively using the DART model. After that, SIF escape probabilities (i.e., 'real' f esc ) were calculated according to Equation (2) in Section 2.2.1. Figure 2 shows the simulated TOC SIF and total emitted SIF, as well as corresponding SIF escape probabilities. It appears that varying input parameters for a given sun direction greatly impacts the simulated TOC SIF, total emitted SIF, and corresponding SIF escape probabilities in all simulations. The variability of canopy structure (i.e., FVC) has a larger impact on TOC SIF (Figure 2a) than total emitted SIF (Figure 2b), which then induced a large variability of SIF escape probabilities from 0.2 to 0.4 ( Figure 2c). Even though leaf chlorophyll content has a considerable impact on TOC SIF and total emitted SIF, respectively, SIF escape probabilities remain relatively stable if leaf chlorophyll content varies from 15 to 60 µg/cm 2 (Figure 2f). On the other hand, with the increase of the solar zenith angle from 20° to 70°, TOC SIF and total emitted SIF decrease at most angles, except from 20° to 30° where the total emitted SIF slightly increases in most cases (Figure 2b,e). At the same time, the corresponding SIF escape probabilities remain relatively stable for all simulations between 30° and 70° (Figure 2c,f).

Simulation of 'Real' f esc in Discontinuous Forest Canopies
Using the variable input parameters listed in Table 2, TOC SIF and total emitted SIF were simulated respectively using the DART model. After that, SIF escape probabilities (i.e., 'real' f esc ) were calculated according to Equation (2) in Section 2.2.1. Figure 2 shows the simulated TOC SIF and total emitted SIF, as well as corresponding SIF escape probabilities. It appears that varying input parameters for a given sun direction greatly impacts the simulated TOC SIF, total emitted SIF, and corresponding SIF escape probabilities in all simulations. The variability of canopy structure (i.e., FVC) has a larger impact on TOC SIF (Figure 2a

Estimation of f esc and SIFtot Using Simulated Canopy Interceptance and fAPARleaf
io and fAPARleaf were first simulated using methods described in Sections 2.2.2 and 2.2.3, and then used to estimate SIF escape probabilities using SIT described in Section 2.1. Figure 3 illustrates the variation of io and fAPARleaf with the increase of FVC, for Cab = 35 µg/cm 2 , background reflectance = Ref-2. Both io and fAPARleaf increase with solar zenith angle. It appeared that io is usually higher than fAPAR for a solar zenith angle larger than 30°.  Figure 4 shows the relative errors of SIF escape probabilities that are estimated using SIT with DART-simulated values (i.e., 'real' values). It appears that FVC generally has a much larger impact on the estimation accuracy of SIF escape probabilities than leaf chlorophyll content. For example, in Figure 4a, varying FVC greatly influences the estimated SIF escape probability with relative errors up to about 180%. In comparison, also in Figure 4a, the relative errors caused by the variation of leaf 3.2. Estimation of f esc and SIF tot Using Simulated Canopy Interceptance and fAPAR leaf i o and fAPAR leaf were first simulated using methods described in Sections 2.2.2 and 2.2.3, and then used to estimate SIF escape probabilities using SIT described in Section 2.1. Figure 3 illustrates the variation of i o and fAPAR leaf with the increase of FVC, for Cab = 35 µg/cm 2 , background reflectance = Ref-2. Both i o and fAPAR leaf increase with solar zenith angle. It appeared that i o is usually higher than fAPAR for a solar zenith angle larger than 30 • . Figure 4 shows the relative errors of SIF escape probabilities that are estimated using SIT with DART-simulated values (i.e., 'real' values). It appears that FVC generally has a much larger impact on the estimation accuracy of SIF escape probabilities than leaf chlorophyll content. For example, in Figure 4a, varying FVC greatly influences the estimated SIF escape probability with relative errors up to about 180%. In comparison, also in Figure 4a, the relative errors caused by the variation of leaf chlorophyll content were no more than 50%. Furthermore, the varying background spectrum has a considerable impact on the estimation of SIF escape probabilities. With a non-reflecting ground, the relative errors between estimation and simulation were no more than 10%. As for the effect of the solar zenith angle, it is clear that the SIF escape probabilities estimated at high solar zenith angles always have smaller relative errors than those estimated at lower solar zenith angles.

Estimation of f and SIFtot Using Simulated Canopy Interceptance and fAPARleaf
io and fAPARleaf were first simulated using methods described in Sections 2.2.2 and 2.2.3, and then used to estimate SIF escape probabilities using SIT described in Section 2.1. Figure 3 illustrates the variation of io and fAPARleaf with the increase of FVC, for Cab = 35 µg/cm 2 , background reflectance = Ref-2. Both io and fAPARleaf increase with solar zenith angle. It appeared that io is usually higher than fAPAR for a solar zenith angle larger than 30°.  Figure 4 shows the relative errors of SIF escape probabilities that are estimated using SIT with DART-simulated values (i.e., 'real' values). It appears that FVC generally has a much larger impact on the estimation accuracy of SIF escape probabilities than leaf chlorophyll content. For example, in Figure 4a, varying FVC greatly influences the estimated SIF escape probability with relative errors up to about 180%. In comparison, also in Figure 4a, the relative errors caused by the variation of leaf chlorophyll content were no more than 50%. Furthermore, the varying background spectrum has a considerable impact on the estimation of SIF escape probabilities. With a non-reflecting ground, the relative errors between estimation and simulation were no more than 10%. As for the effect of the solar zenith angle, it is clear that the SIF escape probabilities estimated at high solar zenith angles always have smaller relative errors than those estimated at lower solar zenith angles. We compared the SIF escape probabilities estimated using canopy interceptance as input and those using fAPARleaf as input (i.e., comparison between Figure 4a and Figure 4c, or Figure 4b and Figure 4d). No apparent differences were observed in most cases. For the estimations using high FVC (i.e., FVC equal to 0.8), however, the estimation accuracy of SIF escape probability using fAPARleaf as input was slightly better than using canopy interceptance as input.
We also compared the estimated SIF escape probabilities using NIRT as input and using NIRv as input that was mentioned in Table 1 (i.e., comparison between Figure 4a and Figure 4b, or Figure 4c and Figure 4d). The accuracy of the SIF escape probability was better when using NIRv than using NIRT. This was especially true with high FVC (i.e., FVC equal to 0.8) as illustrated by the comparison of the estimation accuracy between 0% and 50% with NIRv ( Figure 4b) compared to the accuracy between 50% and 180% with NIRT ( Figure 4a). ) with reference values that were simulated by DART (i.e., ). The calculation formula of was shown in each subplot.
in the first column (a,c) were estimated using NIRT; while in the second column (b,d) using NIRv.
in the first row (a,b) were estimated using io; while in the second row (c,d) using fraction of absorbed photosynthetically active radiation (fAPAR). BS means black soil (i.e., non-reflecting background). Relative error was calculated as ( − )/ .
We further calculated total emitted SIF ( ) using the estimated SIF escape probabilities ( ) mentioned above. Here only the results of °• were used to estimate (corresponding to the in Figure 4b).
in Figure 5a was calculated using specific SIF escape probabilities ( ) that were estimated at varying solar zenith angles from 20° to 70°. While for Figure 5b,c, was calculated using constant SIF escape probabilities ( ) estimated at a fixed solar zenith angle (i.e., SZA = 70° and SZA = 20° for (b) and (c), respectively). In other words, SIF escape probabilities ( ) estimated at SZA = 70° or SZA = 20° were used to assess other cases when solar zenith angles change from 20° to 70°. It is clear that estimated at high solar zenith angles produces higher We compared the SIF escape probabilities estimated using canopy interceptance as input and those using fAPAR leaf as input (i.e., comparison between Figure 4a,c, or Figure 4b,d). No apparent differences were observed in most cases. For the estimations using high FVC (i.e., FVC equal to 0.8), however, the estimation accuracy of SIF escape probability using fAPAR leaf as input was slightly better than using canopy interceptance as input.
We also compared the estimated SIF escape probabilities using NIR T as input and using NIR v as input that was mentioned in Table 1 (i.e., comparison between Figure 4a,b, or Figure 4c,d). The accuracy of the SIF escape probability was better when using NIR v than using NIR T . This was especially true with high FVC (i.e., FVC equal to 0.8) as illustrated by the comparison of the estimation accuracy between 0% and 50% with NIR v (Figure 4b) compared to the accuracy between 50% and 180% with NIR T (Figure 4a).
We further calculated total emitted SIF (SIF SIT tot ) using the estimated SIF escape probabilities ( f esc SIT ) mentioned above. Here only the results of NIR v i• ·ω λ were used to estimate SIF SIT tot (corresponding to the f esc SIT in Figure 4b). SIF SIT tot in Figure 5a was calculated using specific SIF escape probabilities ( f esc SIT ) that were estimated at varying solar zenith angles from 20 • to 70 • . While for Figure 5b

Estimation of f esc and SIFtot Using Remotely Sensed Canopy Interceptance
Canopy interceptance (io) was estimated using the method described in Section 2.3 with leaf area index (LAI) and clumping index (CI), which were then used to estimate SIF escape probability based on SIT. Figure 6 shows the comparison between estimated and simulated results (i.e., 'real' f esc ) of SIF escape probability. Figure 6a shows the estimated results of SIF escape probability using DART-simulated (i.e., 'real') LAI and CI as input for the estimation of canopy interceptance. We found that the estimation accuracy of SIF escape probability was markedly improved (maximum relative error less than 30%) when compared to the estimated results using simulated (i.e., 'real') canopy interceptance (maximum relative error up to 50%, see Figure 4b), except for the cases with a non-reflecting background where the SIF escape probabilities were generally underestimated.
Additional errors can be introduced into the estimation of SIF escape probability when remotely sensed data of LAI and CI were used as input compared to real LAI and CI values. To explore this impact, artificial errors (−50% and +50% constant value error of LAI*CI) were generated and then were used to estimate SIF escape probability (Figure 6b,c). Introducing a −50% error of LAI*CI, SIF escape probabilities were overestimated in most simulations, especially for the cases with high FVC and background mostly covered by vegetation (Ref-2). We noted that the introduction of a +50% error on LAI*CI has no apparent impact on the estimation of SIF escape probability, with only a slight underestimation per simulation.

Estimation of f esc and SIF tot Using Remotely Sensed Canopy Interceptance
Canopy interceptance (i o ) was estimated using the method described in Section 2.3 with leaf area index (LAI) and clumping index (CI), which were then used to estimate SIF escape probability based on SIT. Figure 6 shows the comparison between estimated and simulated results (i.e., 'real' f esc ) of SIF escape probability. Figure 6a shows the estimated results of SIF escape probability using DART-simulated (i.e., 'real') LAI and CI as input for the estimation of canopy interceptance. We found that the estimation accuracy of SIF escape probability was markedly improved (maximum relative error less than 30%) when compared to the estimated results using simulated (i.e., 'real') canopy interceptance (maximum relative error up to 50%, see Figure 4b), except for the cases with a non-reflecting background where the SIF escape probabilities were generally underestimated.
impact, artificial errors (−50% and +50% constant value error of LAI*CI) were generated and then were used to estimate SIF escape probability (Figure 6b,c). Introducing a −50% error of LAI*CI, SIF escape probabilities were overestimated in most simulations, especially for the cases with high FVC and background mostly covered by vegetation (Ref-2). We noted that the introduction of a +50% error on LAI*CI has no apparent impact on the estimation of SIF escape probability, with only a slight underestimation per simulation. Additional errors can be introduced into the estimation of SIF escape probability when remotely sensed data of LAI and CI were used as input compared to real LAI and CI values. To explore this impact, artificial errors (−50% and +50% constant value error of LAI*CI) were generated and then were used to estimate SIF escape probability (Figure 6b,c). Introducing a −50% error of LAI*CI, SIF escape probabilities were overestimated in most simulations, especially for the cases with high FVC and background mostly covered by vegetation (Ref-2). We noted that the introduction of a +50% error on LAI*CI has no apparent impact on the estimation of SIF escape probability, with only a slight underestimation per simulation.
To further evaluate the influence of remotely sensed canopy interceptance on the estimation of SIF escape probability, we calculated SIF SIT tot (Figure 7) using the estimated SIF escape probabilities (i.e., SIF SIT tot = SIF obs / f esc SIT ) mentioned in Figure 6. When compared with the DART-simulated reference values (SIF DART tot ), additional errors of LAI*CI reduce the estimation accuracy of SIF SIT tot (and f esc SIT ) (R 2 reduced from 0.91 to 0.66 and 0.88). Further, underestimated LAI*CI produces a larger influence on the estimation accuracy (R 2 = 0.66) of SIF SIT tot (and f esc SIT ) than overestimated LAI*CI (R 2 = 0.88).
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 16 Figure 6. Comparison of SIF escape probability estimated by SIT and using leaf area index (LAI) and clumping index (CI) as input for the estimation of canopy interceptance (i.e., ) with reference values that were simulated by DART (i.e., ). 'Real' LAI*CI was used to plot (a). Additional −50% and +50% errors of LAI*CI were introduced to plot (b,c). Relative error was calculated as ( − )/ .
To further evaluate the influence of remotely sensed canopy interceptance on the estimation of SIF escape probability, we calculated ( Figure 7) using the estimated SIF escape probabilities (i.e., = / ) mentioned in Figure 6. When compared with the DART-simulated reference values ( ), additional errors of LAI*CI reduce the estimation accuracy of (and ) (R 2 reduced from 0.91 to 0.66 and 0.88). Further, underestimated LAI*CI produces a larger influence on the estimation accuracy (R 2 = 0.66) of (and ) than overestimated LAI*CI (R 2 = 0.88).

Estimation of f esc Using fAPARcanopy and fAPARleaf
SIF escape probabilities were also estimated using two types of fAPAR (i.e., fAPARcanopy = fAPARleaf+woody and fAPARleaf) whose simulation methods were described in Section 2.2.3. Figure 8 shows that the values of SIF escape probability estimated using fAPARcanopy as input were smaller than that using fAPARleaf as input in all simulations. The estimation accuracy of SIF escape probability was improved only for the case with high FVC (FVC equal to 0.8), and the estimation accuracy was

Estimation of f esc Using fAPAR canopy and fAPAR leaf
SIF escape probabilities were also estimated using two types of fAPAR (i.e., fAPAR canopy = fAPAR leaf+woody and fAPAR leaf ) whose simulation methods were described in Section 2.2.3. Figure 8 shows that the values of SIF escape probability estimated using fAPAR canopy as input were smaller than that using fAPAR leaf as input in all simulations. The estimation accuracy of SIF escape probability was improved only for the case with high FVC (FVC equal to 0.8), and the estimation accuracy was decreased for other simulation cases. estimated with additional errors of LAI*CI: −50% error of LAI*CI for (b); +50% error of LAI*CI for (c).

Estimation of f esc Using fAPARcanopy and fAPARleaf
SIF escape probabilities were also estimated using two types of fAPAR (i.e., fAPARcanopy = fAPARleaf+woody and fAPARleaf) whose simulation methods were described in Section 2.2.3. Figure 8 shows that the values of SIF escape probability estimated using fAPARcanopy as input were smaller than that using fAPARleaf as input in all simulations. The estimation accuracy of SIF escape probability was improved only for the case with high FVC (FVC equal to 0.8), and the estimation accuracy was decreased for other simulation cases.

DART-Based Simulation of f esc in Discontinuous Forest Canopies
Although several 3-D SIF models can simulate canopy-leaving SIF [16][17][18], SIF escape probability can be simulated directly only by 3-D models that simulate the SIF emitted per unit area (i.e., per leaf). This is especially important for sparse 3-D vegetation scenes [12]. Since it simulates the total and SIF 3-D radiative budgets, we used the DART model ( Figure S2) to simulate the total SIF emitted by the canopy simultaneously with canopy-leaving SIF radiance in a sparse 3-D vegetation scenes (see methods in Section 2.2 and results in Section 3.1). Furthermore, we also simulated canopy interceptance (i o ) and fAPAR for the same 3-D scenes (Figure 3), which enabled the calculation of escape probabilities using SIT.
Recently, researchers from [34] developed a relatively simple SIF model based on escape and recollision probability (FluorRTER), which is capable to simulate canopy-leaving and total emitted SIF in 3-D discontinuous canopies. Compared to the DART model, however, the scenes simulated by FluorRTER are much simpler. Tree crowns are cylinders, ellipsoids, or conics, and the heterogeneity within tree crowns is not considered, just like its parent model [35]. This missing complexity is important, as woody elements influence forest reflectance and SIF [25,29,36,37]. With the simulation of fAPAR of leaves and woody elements using the DART model, respectively, we evaluated the impact of woody elements on the simulation and estimation of SIF escape probabilities (see Section 3.4). Although the fine structure of the simulation scene is well described in DART, the computational cost of simulating the 3-D radiative budget of very large regions is too large for practical applications. Then, spatial sampling of 3-D radiative budget should be considered in order to keep reasonable computer time and memory according to the complexity of canopy. Here, we did not meet such limitations because we conducted a theoretical analysis of canopy-leaving and total SIF for assessing the SIF escape probability of a complex but small canopy. It is interesting to note that the very recent bi-directional Monte Carlo mode of DART, called DART-Lux, decreases computer time and memory of simulations of bottom and top of the atmosphere SIF radiance images by a hundredfold [38]. This mode opens new perspectives but must be combined with the standard DART mode because it does not simulate 3-D radiative budgets. Moreover, DART-based simulation can also be used to evaluate and validate other models with simplified canopy structure that are more suitable for global application, such as the FluorRTER model.

Influence of Background and Solar Zenith Angle on the Estimation of SIF Escape Probability
We compared the simulation values of SIF escape probability (as 'real' value) with the estimated values using SIT (Figure 4). The estimation accuracy of SIF escape probability is quite high with a non-reflecting background, both with dense forest (FVC = 0.8) and sparse forest (FVC = 0.3), which means that SIT can be directly used in the estimation of SIF escape probability for forest area with a non-reflecting background. However, in the presence of realistic background instead of the non-reflecting background, the estimation accuracy greatly decreased for each simulation setup, especially for the sparse scene (FVC = 0.3) and for a vegetation-covered background . This is because the reflecting background contaminates the NIR v used to estimate the SIF escape probability. This result agrees with previous studies [12,22] and the theoretical limit (i.e., black soil problem) of SIT [19,20,39].
Recently, NIR v (NIR v = NIR T × NDVI) was proposed to replace NIR T for the estimation of SIF escape probability using SIT [12,14], where the use of NDVI is intended to reduce the contamination of soil reflectance on canopy reflectance. Our results confirmed the feasibility of NIR v for specific simulation cases where the background was not totally covered by vegetation. However, the estimation accuracy of SIF escape probability was not improved with the use of NDVI for the simulation case where the background was mostly covered by vegetation (i.e., . Considering the seasonal variation of understory LAI in sparse forests [31,40], the use of NDVI for the estimation of SIF escape probability is still questionable in summer when the background is mostly covered by vegetation (i.e., high understory LAI value).
With regards to the effect of the solar zenith angle on the estimation of SIF escape probability, Figure 4 shows that the accuracies of SIF escape probability estimated at a high solar zenith angle are higher than those estimated at a low solar zenith angle in most studied cases. This phenomenon may be due to the potential link between solar zenith angle and background reflectance. With the solar zenith angle increasing from 20 • to 70 • , more incident photons will be absorbed by the tree crown, and fewer incident photons will be absorbed by the background, which then will weaken the impact of background reflectance on the estimation of SIF escape probability. The above phenomenon was further confirmed by Figure 5b where the estimation accuracy of total emitted SIF was clearly improved when using a constant value of SIF escape probability estimated at 70 • (R 2 increased from 0.85 to 0.95). Conversely, the estimation accuracy of total emitted SIF slightly decreased (R 2 decreased from 0.85 to 0.81) when using a constant value of SIF escape probability estimated at 20 • (Figure 5c). It is worth noting that these findings were based on the assumption that other input parameters remain constant when the solar zenith angle changes from 20 • to 70 • (i.e., over the season). In reality, it is not the case, because canopy structure and leaf optical properties change over the season [1,31,41]. Fortunately, Figure 2 shows that the impact of leaf optical properties (represented by Cab) on SIF escape probability is slight in the far-red region, which agrees with results obtained with the SCOPE model ( Figure 4 in [22]).

Influence of Remotely Sensed i o and fAPAR on the Estimation SIF Escape Probability
In general, remotely sensed canopy interceptance (Figure 6a) led to a relatively accurate SIF escape probability compared to the simulated canopy interceptance (i.e., 'real' value, see Figure 4b).
Interestingly, the use of remotely sensed canopy interceptance slightly improved the estimation accuracies for configurations with the natural background spectrum. This phenomenon was explained by the joint influences of canopy interceptance and canopy reflectance on the estimation of SIF escape probability, which we term the contamination compensation effect. Firstly, SIF escape probability was generally overestimated by the use of canopy reflectance (NIR v or NIR T ), because canopy reflectance is contaminated by the reflectance of the background (see Figure 4). Then, the use of overestimated canopy interceptance ( Figure S3a-c) partly reduces the estimation errors caused by the contaminated canopy reflectance, especially for the sparse forest canopies. On the contrary, for these configurations with a non-reflecting background spectrum, because there is no more contamination caused by the background reflectance, the estimation accuracy of SIF escape probability was mostly reduced (Figure 6a). However, it should be noted that the above analysis was based on the assumption that G(θ) was fixed to 0.5, as for global-scale remote sensing applications [13,14,25,33], and LAI and CI were set to simulation values (i.e., 'real' value).
The introduction of artificial constant value errors on LAI*CI has a remarkable influence on the estimation accuracy of SIF escape probability (and total emitted SIF), with a clear overestimate of SIF escape probabilities if LAI*CI is underestimated (Figure 6b). To better explain this trend, canopy interceptance was plotted using LAI*CI with corresponding artificial errors ( Figure S3a-c). It demonstrated that an underestimated LAI*CI implies an underestimated canopy interceptance, which amplifies the overestimation induced by the contaminated canopy reflectance. An overestimated LAI*CI tends to produce an overestimated canopy interceptance, which is somewhat compensated by the contaminated canopy reflectance.
The errors of satellite-retrieved LAI and CI products, including theoretical and physical errors, are inevitable [42]. These will be more significant for forest biomes where the accurate estimation of LAI remains challenging with an uncertainty relatively higher than for other biomes [43][44][45]. In a recent study, Zhang et al. [14] demonstrated that total emitted SIF, estimated using satellite-retrieved LAI and CI products, was better correlated with GPP for all C3 sites combined, but was less correlated for forest sites. The above phenomenon could be partly attributed to the high error and uncertainty of satellite-retrieved LAI and CI products. Our result (see Figure 7) also demonstrated the high uncertainty of estimation of canopy interceptance in the forest using LAI and CI, especially for underestimated LAI*CI. Furthermore, G(θ) was fixed to 0.5 for all biomes, which can also be a potential source of error because it only holds for homogeneous canopy with spherical leaf orientation distribution [25,33]. Overall, the estimation of canopy interceptance, and the SIF escape probability or total emitted SIF, using satellite-retrieved LAI and CI products is still questionable and needs further investigation.
Owing to the high uncertainty of estimated canopy interceptance using satellite-retrieved LAI and CI products, Zeng et al. [12] proposed to use fAPAR to approximate canopy interceptance instead of canopy interceptance. Our simulation results support this approach since it leads to comparable and even better estimation of SIF escape probability than the method using canopy interceptance ( Figure 3). However, fAPAR for whole canopy (i.e., fAPAR canopy ) generally can be divided into fAPAR for photosynthetic elements (i.e., fAPAR leaf ) and fAPAR for non-photosynthetic elements (i.e., fAPAR woody ), respectively [46,47]. The impact of fAPAR woody is neglected in most satellite-retrieved fAPAR products [48]. As a result, the direct use of satellite-retrieved fAPAR products (mostly fAPAR canopy ) to estimate SIF escape probability still needs further investigation [12,49]. Our simulation results illustrated that using fAPAR canopy will induce an underestimation of SIF escape probability when compared to fAPAR leaf (Figure 8). This underestimation can be a potential advantage for the sparse forest where the SIF escape probabilities are generally overestimated due to the contamination of background, providing an additional mechanism behind the contamination compensation effect. Still, more measurement data from stand to global scale are needed to explore this finding and develop a robust method to separate the various canopy layers.

Conclusions
This study evaluated the feasibility of SIF escape probability estimation methods in discontinuous forest canopies. In general, these estimation methods based on spectral invariants theory have good accuracy for relatively dense forest canopies. The use of NDVI can partly improve the low estimation accuracy for sparse forest canopies. However, if the forest background is partly or totally covered by vegetation, NDVI failed to eliminate the effects of background reflectance on the estimation of SIF escape probability, because NDVI failed to totally reduce the contamination of background reflectance on canopy reflectance. SIF escape probability estimated at a high solar zenith angle will have a higher accuracy than those estimated at a low solar zenith angle in most studied cases due to the potential link between solar zenith angle and background reflectance. For global-scale remote sensing applications, spectral invariant canopy interceptance can be remotely sensed using the current satellite LAI and CI products. However, additional errors will be introduced with the use of satellite products, especially when the product of LAI and CI was underestimated. Our simulation results also show that canopy interceptance can be replaced by fAPAR with a comparable estimation accuracy, and the direct use of fAPAR canopy may have a potential advantage than fAPAR leaf in sparse forest canopies due to a contamination compensation mechanism. To summarize, our results improved our understanding of SIF escape probability estimation methods in discontinuous forest canopies. However, future work is required to evaluate these findings using measurement data. Such work could explore OCO-2 or TROPOMI SIF retrievals or field observations. Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/23/3962/s1, Figure S1: Three background reflectance spectra used for DART simulations (Ref-1: bare soil; Ref-2: soil is partly covered by vegetation; Ref-3: soil is totally covered by vegetation), Figure S2: DART-simulated radiative budgets (Intercepted, Scattered, Absorbed, and Emitted radiance) at each layer of the forest crown (a), emitted SIF at each layer under different simulation settings (b). SIF was simulated at 740 nm. The solar zenith angle was set to 30 • for all simulations, Figure S3: Simulated canopy interceptance (i.e., first row) and corresponding SIF escape probabilities (i.e., second row) using 'real' LAI*CI, −50% error LAI*CI, and +50% error LAI*CI.