Evaluation of the MODIS LAI/FPAR Algorithm Based on 3D-RTM Simulations: A Case Study of Grassland

Uncertainty assessment of the moderate resolution imaging spectroradiometer (MODIS) leaf area index (LAI) and the fraction of photosynthetically active radiation absorbed by vegetation (FPAR) retrieval algorithm can provide a scientific basis for the usage and improvement of this widely-used product. Previous evaluations generally depended on the intercomparison with other datasets as well as direct validation using ground measurements, which mix the uncertainties from the model, inputs, and assessment method. In this study, we adopted the evaluation method based on three-dimensional radiative transfer model (3D RTM) simulations, which helps to separate model uncertainty and other factors. We used the well-validated 3D RTM LESS (large-scale remote sensing data and image simulation framework) for a grassland scene simulation and calculated bidirectional reflectance factors (BRFs) as inputs for the LAI/FPAR retrieval. The dependency between LAI/FPAR truth and model estimation serves as the algorithm uncertainty indicator. This paper analyzed the LAI/FPAR uncertainty caused by inherent model uncertainty, input uncertainty (BRF and biome classification), clumping effect, and scale dependency. We found that the uncertainties of different algorithm paths vary greatly (−6.61% and +84.85% bias for main and backup algorithm, respectively) and the “hotspot” geometry results in greatest retrieval uncertainty. For the input uncertainty, the BRF of the near-infrared (NIR) band has greater impacts than that of the red band, and the biome misclassification also leads to nonnegligible LAI/FPAR bias. Moreover, the clumping effect leads to a significant LAI underestimation (−0.846 and −0.525 LAI difference for two clumping types), but the scale dependency (pixel size ranges from 100 m to 1000 m) has little impact on LAI/FPAR uncertainty. Overall, this study provides a new perspective on the evaluation of LAI/FPAR retrieval algorithms.


Introduction
Leaf area index (LAI), defined as half of the total green leaf area of per unit horizontal ground area, is a basic parameter for measuring the vegetation canopies [1,2]. This variable plays a key roles in hydrology, biogeochemistry, and ecosystem models that connect vegetation to the climate observing system through the carbon, water cycles, and radiation [3]. Fraction of photosynthetically

MODIS LAI/FPAR Retrieval Algorithm
The MODIS LAI/FPAR retrieval algorithm consists of a main algorithm based on the radiative transfer equation (RTE) and a backup algorithm using the relationship between vegetation index and LAI/FPAR. The retrieval algorithm exploits the spectral information content of MODIS surface reflectances at up to 7 spectral bands (band 1: 620-670 nm; band 2: 841-876 nm; band 3: 459-479 nm; band 4: 545-565 nm; band 5: 1230-1250 nm; band 6: 1628-1652 nm; band 7: 2105-2155 nm) [4,8]. Inputs of this algorithm include BRFs at red and near-infrared (NIR) bands (band 1 and 2), their uncertainties, sun-sensor geometry (SZA: solar zenith angle, SAA: solar azimuth angle, VZA: view zenith angle, VAA: view azimuth angle), and a biome classification map. Note that in the current algorithm version, different biome types use different RT models. Herbaceous biomes (B1: grasses and cereal crops; B2: shrubs; B3: broadleaf crops;) were modelled using 1D RT due to the good continuity of the grass distribution and in consideration of the computational efficiency. Savannas (B4) were modelled by a stationary Poisson germ-grain stochastic process (so called stochastic radiative transfer (SRT) model) [47,48]. Forest biomes (B5: evergreen broadleaf forests; B6: deciduous broadleaf forests; B7: evergreen needleleaf forests; and B8: deciduous needleleaf forests) were based on a 3D RTM (3D structures were represented by columns uniformly (deterministically) spaced on the ground). With these RTMs, the science team constructed an LAI/FPAR main algorithm based on angular information, biome type, and spectral information in which the mean and standard deviation values of the LAI and FPAR selected in the spectral retrieval space are reported for retrieval value and its uncertainty. The main look up table (LUT)-based algorithm was designed as follows. Firstly, the main algorithm evaluates a weight coefficient as a function of sun-sensor geometry, wavelength, and LAI by using a field-tested canopy reflectance model. Then it calculates the BRFs by using the weight coefficient and the same model [4,8]. The algorithm tests the eligibility of a canopy radiation model to generate the LUT file where a subset of coefficients is satisfied within a given accuracy [9]. The given atmosphere-corrected BRFs are then compared with the modeled BRFs, which are stored in the biome-specific LUT files. Finally, all candidates of LAI/FPAR are used to calculate the mean values and uncertainty of the retrieval [9]. In the case of highly dense canopies, reflectance will be saturated and insensitive to changes in canopy properties. Therefore, LAI and FPAR values acquired under saturated conditions are less reliable than those generated by unsaturated BRFs. When the main algorithm fails to localize a solution, the backup algorithm is used to retrieve values through an empirical relationship between the normalized difference vegetation index (NDVI) and the canopy LAI/FPAR [11,21]. Such retrievals are flagged in the algorithm path quality assessment (QA) variable [8], which consists of two values for the main algorithm and two values for the backup algorithm (from high quality to low): the main algorithm without saturation (QA = 0), the main algorithm with saturation (QA = 1), the backup algorithm due to sun-sensor geometry (QA = 3), and the backup algorithm due to other reasons (QA = 4) [9,14,49].

Three-Dimensional Grassland Scene Simulation
We used the newly proposed but well validated 3D RT model LESS to simulate the interaction between the solar radiation and landscape elements based on the spectral response functions (SRFs) (from ENVI software) of MODIS and calculated the scene BRFs [35,45]. LESS simulates BRFs by a weighted forward photon tracing method as well as simulated energy transfer and generates images by a backward path tracing method [35]. Qi et al. [35] described the comparison between BRFs simulated by LESS and average BRF results from other models (e.g., SPRINT3, RAYTRAN, and RAYSPREAD) over several different homogeneous and heterogeneous canopies from the RAMI website to evaluate the accuracy of LESS.
The input parameters of LESS include 3D landscape elements, optical properties, and sun-sensor geometries. The simulated scenes are covered by grass (Johnson grass) and its component spectra were obtained from the LOPEX93 dataset on the OPTICLEAF website [50]. The soil (grayish brown loam) spectra were selected from the soil spectral library in ENVI software, and the transmittance of the soil is 0 ( Figure 1). Then we calculated the two MODIS bands (red: band 1 and NIR: band 2) reflectance and transmittance (see Table 1) using SRFs (the shaded part of Figure 1) of the MODIS sensor by the following equation [51]: where, R and T are MODIS band reflectance and transmittance, respectively. The R λ and T λ are mean narrow-band reflectance and transmittance derived from the spectral curves. The S λ is the SRF value of the MODIS sensor. λ is the value of wavelength, which has a specific upper (λ max : red = 670 µm, NIR = 876 µm) and lower (λ min : red = 620 µm, NIR = 841 µm) limit for each band.

Three-dimensional Grassland Scene Simulation
We used the newly proposed but well validated 3D RT model LESS to simulate the interaction between the solar radiation and landscape elements based on the spectral response functions (SRFs) (from ENVI software) of MODIS and calculated the scene BRFs [35,45]. LESS simulates BRFs by a weighted forward photon tracing method as well as simulated energy transfer and generates images by a backward path tracing method [35]. Qi et al. [35] described the comparison between BRFs simulated by LESS and average BRF results from other models (e.g., SPRINT3, RAYTRAN, and RAYSPREAD) over several different homogeneous and heterogeneous canopies from the RAMI website to evaluate the accuracy of LESS.
The input parameters of LESS include 3D landscape elements, optical properties, and sunsensor geometries. The simulated scenes are covered by grass (Johnson grass) and its component spectra were obtained from the LOPEX93 dataset on the OPTICLEAF website [50]. The soil (grayish brown loam) spectra were selected from the soil spectral library in ENVI software, and the transmittance of the soil is 0 ( Figure 1). Then we calculated the two MODIS bands (red: band 1 and NIR: band 2) reflectance and transmittance (see Table 1) using SRFs (the shaded part of Figure 1) of the MODIS sensor by the following equation [51]: where, R and T are MODIS band reflectance and transmittance, respectively. The and are mean narrow-band reflectance and transmittance derived from the spectral curves. The is the SRF value of the MODIS sensor.
is the value of wavelength, which has a specific upper ( : red = 670 , NIR = 876 ) and lower ( : red = 620 , NIR = 841 ) limit for each band.   The 3D landscape elements were created with the third-party software OnyxTree, which uses the calculated reflectance and transmittance (Table 1) to make a grass 3D model (obj format file). As shown in Figure 2, we created nine randomly distributed grasslands with different LAIs (0.25, 0.50, 0.75, 1.0, 1.25, 1.5, 2.5, 3.5, and 4.5) using the LESS and grass 3D model. Moreover, LESS calculates FPAR by performing a band integration of the PAR between 380 nm and 710 nm and dividing by the incident radiation (slightly different from MODIS for which the wavelength interval is 400-800 nm) based on the LESS simulation of the collision of photons and the transfer of energy. In addition, to match the canopy structure of grasses in the MODIS LAI/FPAR retrieval algorithm (all organs other than leaves are ignored), only foliage is present in the scene. There is also only direct radiation in these scenes. The size of these scenes is 500 m × 500 m, which matches the spatial resolution of the MODIS LAI/FPAR products.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 18 represent the SRFs of the MODIS sensor in the red (620-670 ) and NIR (841-876 ) bands, respectively. Table 1. Broad-band reflectance and transmittance of grass and soil used in this study. R and T are abbreviations for broad-band reflectance and transmittance, respectively.

R (red) T (red) R (NIR) T (NIR)
Johnson grass 0.0738 0.0577 0.4276 0.4607 Grayish brown loam 0.1755 0 0.3021 0 The 3D landscape elements were created with the third-party software OnyxTree, which uses the calculated reflectance and transmittance (Table 1) to make a grass 3D model (obj format file). As shown in Figure 2, we created nine randomly distributed grasslands with different LAIs (0.25, 0.50, 0.75, 1.0, 1.25, 1.5, 2.5, 3.5, and 4.5) using the LESS and grass 3D model. Moreover, LESS calculates FPAR by performing a band integration of the PAR between 380 nm and 710 nm and dividing by the incident radiation (slightly different from MODIS for which the wavelength interval is 400-800 nm) based on the LESS simulation of the collision of photons and the transfer of energy. In addition, to match the canopy structure of grasses in the MODIS LAI/FPAR retrieval algorithm (all organs other than leaves are ignored), only foliage is present in the scene. There is also only direct radiation in these scenes. The size of these scenes is 500 m × 500 m, which matches the spatial resolution of the MODIS LAI/FPAR products.

Experimental Design
We utilized the standard deviation of all LAI/FPAR candidates (StdLAI and StdFPAR), the retrieval index (RI), and the relative and absolute LAI/FPAR differences as the indicators of LAI/FPAR uncertainty. According to the uncertainty theory, StdLAI and StdFPAR are the standard deviations of all acceptable LAI/FPAR solutions in the LUT, which are the function of both the input uncertainty (biome type and BRF uncertainty) and model uncertainty [4,8]. StdLAI and StdFPAR have been proven and evaluated as quality metrics for MODIS LAI/FPAR products [28,52]. However, these two metrics have limitations due to the regularization introduced by the LUT algorithm and are artificially lowered at large LAIs [8]. Therefore, in this paper, we have also selected the RI (see Equation (2)) as an uncertainty metric, which is defined as the percentage of pixels for which the main RTE-based algorithm generates retrieval results. We note that the RI is used to characterize the overall uncertainty of all pixels [21,25,53], while StdLAI and StdFPAR are used to characterize individual pixel uncertainty.

RI =
Number o f pixels retrieved by the main algorithm Total number o f processed pixels (2) To evaluate the consistency between true LAI/FPAR and MODIS retrievals, the difference between the simulation results of LAI (input to the LESS)/FPAR (output from the LESS) and the LAI/FPAR retrieved by the MODIS retrieval algorithm were used. The relative difference (RD, see Equation (3)) and absolute difference (AD, see Equation (4)), were utilized to quantify any differences.
Based on the uncertainty theory, the retrieval uncertainty is a function of both model and input uncertainty and is embedded in the MODIS algorithm. In this study, we explored the relationship between retrieval uncertainty and the retrieval space, sun-sensor geometry, surface reflectance uncertainty, and biome type uncertainty using the variable-controlling approach (see Table 2). We analyzed the inherent model uncertainty in two steps: 1) analysis of the retrieval space; 2) uncertainty changed with sun-sensor geometry. We obtained 4000 red-NIR BRF pairs by adding normally distributed errors (errors with 5% and 15% standard deviation) to the LESS simulated red and NIR band BRFs (1000: red without uncertainties and NIR with 5% standard deviation, 1000: red without uncertainties and NIR with 15% standard deviation, 1000: NIR without uncertainties and red with 5% standard deviation, and 1000: NIR without uncertainties and red with 15% standard deviation). Then we analyzed the LAI/FPAR uncertainty caused by BRF uncertainty within the 4 groups of samples. In addition, we analyzed uncertainties due to biome type misclassification, which is one of the main factors affecting the LAI/FPAR retrieval accuracy [4,54]. Each red-NIR BRF pair was sequentially combined with each biome type as the inputs for the MODIS LAI/FPAR algorithm. In this experiment, only B1 (grasses and cereal crops) was correct while the remaining seven combinations represented the biome type misclassification cases. Finally, we analyzed the influence of scale dependency and clumping effect ("tree groups") [55] on the uncertainty of LAI/FPAR retrievals. We simulated a randomly distributed 1 km × 1 km scene (Figure 9a-1) and two clumping 1 km × 1 km scenes. Clumping type 1 (CT1, Figure 9a-2) had random clumping and Clumping type 2 (CT2, Figure 9a-3) was half bare ground and half grass. The LAI of the three scenes remained constant and these scenes were downscaled into four 500 m × 500 m scenes, sixteen 250 m × 250 m scenes, and one hundred 100 m × 100 m scenes for the discussion of scale dependency and clumping effect.

Inherent Model Uncertainty
To evaluate the inherent model uncertainty of the MODIS LAI/FPAR retrieval algorithm, we analyzed the effect of the retrieval space and sun-sensor uncertainty, separately. In performing the evaluation of the retrieval space, we paid more attention to the changes in the uncertainty of algorithm paths. While the difference between the LAI/FPAR retrieval and LESS simulations were analyzed when evaluating of the sun-sensor geometry. Figure 3 indicates the variation of LAI/FPAR and its uncertainty in the retrieval space. As we can see, LAI/FPAR is nonlinearly related to surface reflectance (Figure 3a,b), and FPAR is also nonlinearly related to LAI. Moreover, the relationship between LAI/FPAR and its uncertainty (StdLAI and StdFPAR) is also nonlinear. The StdLAI and StdFPAR are very low for lower LAI/FPAR and then increase to the highest values, and then steadily decrease (from the bottom right to the top left of Figure 3d,e) as the LAI/FPAR gets progressively larger (from the bottom right to the top left of Figure 3a,b). It is also obvious that there is a clear division between the saturated (QA = 1) and unsaturated (QA = 0) parts where the LAI/FPAR values are higher in the saturated part (Figure 3c). Compared to the unsaturated part, the bias of LAI (+4.64) and FPAR (+0.631) are high, but the bias of StdLAI (−0.052) and StdFPAR (−0.169) are low in the saturated part. Figure 3d,e show that StdLAI and StdFPAR are relatively small at the boundaries of the area retrieved by the main algorithm due to the regularization of the algorithm [4,8].

Retrieval Uncertainty as a Function of Sun-Sensor Geometry
The relationship between LAI/FPAR uncertainty and sun-sensor geometry are presented in Figures 4 and 5. In the high LAI scene (Figure 4a, LAI = 3.50), the retrieval results of LAI/FPAR show low consistency with the truth (it yields to an overall uncertainty of 20.01% for RD of LAI and 13.96% for RD of FPAR). The main algorithm shows an averaged 6.61% underestimation of LAI, while the backup algorithm results in an averaged 84.85% overestimation of LAI. In this scene, the backup algorithm appears at the "hotspot" geometry and where the difference between SAA and VAA is large. It can also be seen that the large VZA will lead to saturation. Nevertheless, for the low LAI scene (LAI = 0.50), the retrieved LAI/FPAR showed a significant overestimation (+111.86% RD for LAI, +162.50% RD for

Retrieval Uncertainty as a Function of Sun-Sensor Geometry
The relationship between LAI/FPAR uncertainty and sun-sensor geometry are presented in Figure 4 and Figure 5. In the high LAI scene (Figure 4a, LAI = 3.50), the retrieval results of LAI/FPAR show low consistency with the truth (it yields to an overall uncertainty of 20.01% for RD of LAI and 13.96% for RD of FPAR). The main algorithm shows an averaged 6.61% underestimation of LAI, while the backup algorithm results in an averaged 84.85% overestimation of LAI. In this scene, the backup algorithm appears at the "hotspot" geometry and where the difference between SAA and VAA is large. It can also be seen that the large VZA will lead to saturation. Nevertheless, for the low LAI scene (LAI = 0.50), the retrieved LAI/FPAR showed a significant overestimation (+111.86% RD for LAI, +162.50% RD for FPAR) and large uncertainty (StdLAI = 0.285, and StdFPAR = 0.238). Figure 5 shows the same analysis as above, but we controlled the view position and varied the sun position. Comparing Figures 4 and 5, the distribution of LAI and its uncertainty (Figure 4a

Retrieval Uncertainty as a Function of Sun-Sensor Geometry
The relationship between LAI/FPAR uncertainty and sun-sensor geometry are presented in Figure 4 and Figure 5. In the high LAI scene (Figure 4a, LAI = 3.50), the retrieval results of LAI/FPAR show low consistency with the truth (it yields to an overall uncertainty of 20.01% for RD of LAI and 13.96% for RD of FPAR). The main algorithm shows an averaged 6.61% underestimation of LAI, while the backup algorithm results in an averaged 84.85% overestimation of LAI. In this scene, the backup algorithm appears at the "hotspot" geometry and where the difference between SAA and VAA is large. It can also be seen that the large VZA will lead to saturation. Nevertheless, for the low LAI scene (LAI = 0.50), the retrieved LAI/FPAR showed a significant overestimation (+111.86% RD for LAI, +162.50% RD for FPAR) and large uncertainty (StdLAI = 0.285, and StdFPAR = 0.238). Figure 5 shows the same analysis as above, but we controlled the view position and varied the sun position. Comparing Figures 4 and 5, the distribution of LAI and its uncertainty (Figure 4a

Input BRF Uncertainty
Here, we calculated the effects of input BRF uncertainty on the LAI/FPAR retrieval. Figure 6 shows that the uncertainty in the LAI/FPAR in the shadow area of the 15% BRF uncertainty is much larger, which means that larger BRF uncertainty will result in larger LAI/FPAR uncertainty. The StdLAI and StdFPAR due to a 5% BRF uncertainty are close to the StdLAI and StdFPAR due to a 15% BRF in the red band. This is because both 5% and 15% BRF uncertainty in the red band will trigger the backup algorithm with no StdLAI and StdFPAR. Comparing the shadow area in panels (a) and (b), we found that the same level of uncertainty in the NIR band BRF has a greater impact on the retrieval than the red band BRF. The main algorithm was not used in the hotspot (VZA = 0) geometry leading to the absence of both StdLAI and StdFPAR in panel (a).

Input BRF Uncertainty
Here, we calculated the effects of input BRF uncertainty on the LAI/FPAR retrieval. Figure 6 shows that the uncertainty in the LAI/FPAR in the shadow area of the 15% BRF uncertainty is much larger, which means that larger BRF uncertainty will result in larger LAI/FPAR uncertainty. The StdLAI and StdFPAR due to a 5% BRF uncertainty are close to the StdLAI and StdFPAR due to a 15% BRF in the red band. This is because both 5% and 15% BRF uncertainty in the red band will trigger the backup algorithm with no StdLAI and StdFPAR. Comparing the shadow area in panels (a) and (b), we found that the same level of uncertainty in the NIR band BRF has a greater impact on the retrieval than the red band BRF. The main algorithm was not used in the hotspot (VZA = 0) geometry leading to the absence of both StdLAI and StdFPAR in panel (a).

Input Biome Type Uncertainty
Different biome types have different canopy structures, and the MODIS retrieval algorithm uses photon transport theory and the corresponding RT model for different biome types to parameterize the canopy structures (e.g., reflectance and transmittance of leaves, crown shadowing), which form the LUTs of the MODIS retrieval algorithm. To check the sensitivity of the algorithm to biome type, we modified the input biome type for the retrieval algorithm from the correct type (B1: grasses and cereal crops) to incorrect types (B2: shrubs; B3: broadleaf crops; B4: savannas; B5: evergreen broadleaf forests; B6: deciduous broadleaf forests; B7: evergreen needleleaf forests; and B8: deciduous needleleaf forests). As seen from Figures 7 and 8, the retrieval uncertainty is similar when the input biome types are non-forest biomes (B1-B4) with a greater than 59.5% RI for all four biome types except for B2 in the LAI = 4.5 scene. However, the RI gets much lower when the grassland pixel is misclassified into forest biomes. As shown in Figures 7c and 8c, the uncertainties of the retrieved LAI are high at the scene with high LAI (e.g., LAI = 3.5, 4.5). A significant overestimation (+0.727, +1.434 for AD of LAI) in B2, and a significant underestimation (−1.608, −2.344 for AD of LAI) in B3 is also evident. For B5, the RI is high (>69%) but AD of LAI (>1.656) is also high when LAI is relatively high (e.g., LAI = 2.5, 3.5, 4.5). As shown in Figure 8c, the FPAR calculated from the MODIS algorithm is significantly overestimated for all cases except for B4 high LAI scenes, which appear to be underestimated.

Input Biome Type Uncertainty
Different biome types have different canopy structures, and the MODIS retrieval algorithm uses photon transport theory and the corresponding RT model for different biome types to parameterize the canopy structures (e.g., reflectance and transmittance of leaves, crown shadowing), which form the LUTs of the MODIS retrieval algorithm. To check the sensitivity of the algorithm to biome type, we modified the input biome type for the retrieval algorithm from the correct type (B1: grasses and cereal crops) to incorrect types (B2: shrubs; B3: broadleaf crops; B4: savannas; B5: evergreen broadleaf forests; B6: deciduous broadleaf forests; B7: evergreen needleleaf forests; and B8: deciduous needleleaf forests). As seen from Figures 7 and 8, the retrieval uncertainty is similar when the input biome types are non-forest biomes (B1-B4) with a greater than 59.5% RI for all four biome types except for B2 in the LAI = 4.5 scene. However, the RI gets much lower when the grassland pixel is misclassified into forest biomes. As shown in Figure 7c or Figure 8c, the uncertainties of the retrieved LAI are high at the scene with high LAI (e.g., LAI = 3.5, 4.5). A significant overestimation (+0.727, +1.434 for AD of LAI) in B2, and a significant underestimation (−1.608, −2.344 for AD of LAI) in B3 is also evident. For B5, the RI is high (>69%) but AD of LAI (>1.656) is also high when LAI is relatively high (e.g., LAI = 2.5, 3.5, 4.5). As shown in Figure 8c, the FPAR calculated from the MODIS algorithm is significantly overestimated for all cases except for B4 high LAI scenes, which appear to be underestimated.

Impact of Clumping Effect and Scale Dependency
The model scale dependency and clumping effect have attracted much attention from the community in the development of quantitative remote sensing. In this experiment, the model scale dependency refers to the discrepancy between LAI/FPAR uncertainties that are derived from the same algorithm but at different spatial resolutions. The scale dependency determines the adaptive capacity of an algorithm for different pixel size. The clumping effect refers to the discrepancy between retrieved LAI/FPARs with same LAI/FPAR truth but different vegetation spatial distributions. The model nonlinearly and surface heterogeneity together result in the well-known phenomenon called "Inversion first and aggregation later is different from aggregation first and inversion later" [11].
Comparing the algorithm performance at different scales, we found that the MODIS algorithm is nearly scale-invariant from 100 m to 1000 m. Both LAI/FPAR and their uncertainty nearly remain unchanged with increasing pixel size except for the CT2, which shows that the StdLAI and StdFPAR are lower than other scales ( Table 3). The retrieved LAIs for all scenes are less than the LAI truth at 1000 m scale, and the underestimations for Uniform, CT1, and CT2 vegetation distributions are −0.005, −0.846, and −0.525, respectively. Comparing the three clumping scenes, we found that the LAI of a uniform scene is very close to the LAI truth (Figure 9b-1). CT1 shows a significant underestimation, while CT2 shows a significant overestimation except in the 1000 m scale. For FPAR, there is a significant overestimation in all three scenes (Figure 9b-2). At the same spatial resolution, the RI of CT1 is the highest, followed by Uniform, and the lowest is CT2 (Table 3). While the values of StdLAI and StdFPAR are as follows (from small to large): Uniform, CT1, and CT2. The standard deviations of StdLAI and StdFPAR also get larger in this order.

Impact of Clumping Effect and Scale Dependency
The model scale dependency and clumping effect have attracted much attention from the community in the development of quantitative remote sensing. In this experiment, the model scale dependency refers to the discrepancy between LAI/FPAR uncertainties that are derived from the

Impact of Clumping Effect and Scale Dependency
The model scale dependency and clumping effect have attracted much attention from the community in the development of quantitative remote sensing. In this experiment, the model scale dependency refers to the discrepancy between LAI/FPAR uncertainties that are derived from the  same algorithm but at different spatial resolutions. The scale dependency determines the adaptive capacity of an algorithm for different pixel size. The clumping effect refers to the discrepancy between retrieved LAI/FPARs with same LAI/FPAR truth but different vegetation spatial distributions. The model nonlinearly and surface heterogeneity together result in the well-known phenomenon called "Inversion first and aggregation later is different from aggregation first and inversion later" [11]. Comparing the algorithm performance at different scales, we found that the MODIS algorithm is nearly scale-invariant from 100 m to 1000 m. Both LAI/FPAR and their uncertainty nearly remain unchanged with increasing pixel size except for the CT2, which shows that the StdLAI and StdFPAR are lower than other scales ( Table 3). The retrieved LAIs for all scenes are less than the LAI truth at 1000 m scale, and the underestimations for Uniform, CT1, and CT2 vegetation distributions are -0.005, −0.846, and −0.525, respectively. Comparing the three clumping scenes, we found that the LAI of a uniform scene is very close to the LAI truth (Figure 9b-1). CT1 shows a significant underestimation, while CT2 shows a significant overestimation except in the 1000 m scale. For FPAR, there is a significant overestimation in all three scenes (Figure 9b-2). At the same spatial resolution, the RI of CT1 is the highest, followed by Uniform, and the lowest is CT2 (Table 3). While the values of StdLAI and StdFPAR are as follows (from small to large): Uniform, CT1, and CT2. The standard deviations of StdLAI and StdFPAR also get larger in this order.

Discussion
Because of the different sensitivities of LAI/FPAR to surface reflectances, we note that there would be a gap of uncertainty between the saturated part and the unsaturated part [4,8]. However, Figure 3 indicates that for large LAI/FPAR, their theoretical uncertainty is artificially reduced by the method of regularization, which causes the retrieval to have varying degrees of confidence and leads to a problematic evaluation of high LAI/FPAR scenes using the provided StdLAI and StdFPAR. This also places new requirements on future algorithm refinement that the LUT algorithm should be consistent in the saturated case as in the unsaturated case. The LAI/FPAR values estimated by the backup algorithm and calculated by the main algorithm also show significant discontinuity [4,[8][9][10] (Figure 4a). Based on this, we point out that future algorithm refinement should increase the coverage of the main algorithm usage, which will greatly improve the overall accuracy of the product. In addition, according to the 3D RT model, the hotspot means that the radiation field tends to peak around the retro-illumination direction. The results of this study indicate that the uncertainty of the MODIS algorithm in the hotspots is quite large (Figures 4-6), due to which the science team decided not to include additional hotspot parameters since their inclusion would make algorithm calibration difficult [56,57]. We note that this will not cause large problems in the MODIS LAI/FPAR production because of the fact that observations near the hotspot are rare for MODIS. However, this points out a new refinement direction of this algorithm to improve the accuracy of hotspot modeling for other sensors.
As is known, the uncertainty of the inputting BRFs has some influence on the uncertainty of the retrieval algorithm. In particular, our results show that the uncertainty of NIR BRFs has a larger effect on LAI/FPAR uncertainty compared to the red BRFs ( Figure 6). We know that insufficient input information will lead to the "ill-posed" retrieval problem [11]; however, the inputting BRFs of the MODIS operational algorithm are currently only for the red and NIR bands. Therefore, in the future we may try to make use of BRFs in other bands to improve the retrieval accuracy. The MODIS algorithm depends on a priori information about the land surface given by biome type representing the pattern of the architecture of vegetation, as well as patterns of spectral reflectance and transmittance of vegetation [8]. Figures 7 and 8 confirm that the misclassification of biome types with similar structures will result in smaller LAI/FPAR uncertainty, and vice versa [11,58]. This means that the improvement of biome classification accuracy is an efficient way to improve the LAI/FPAR products. Moreover, different biome types also lead to different clumping types. As Figure 9 shows, the underestimation of LAI is significant for two clumping scenes at 1000 m scale. For the other three scales, however, the overestimation of CT2 is due to the backup algorithm retrievals. As our results show, the algorithm only considers the clumping effect at one scale (e.g., B1 is minimal leaf clumping) [4], which can result in large differences in the retrievals; therefore, we suggest that future algorithms consider the clumping effect at more scales (e.g., leaf, branch, and crown).
We note that there are some problems with the way we use LESS to simulate specific scenes and evaluate the MODIS algorithm. First, according to the algorithm, the retrieved LAI/FPAR is a weighted average of the probability values within the error range. Therefore, the probability distribution of LAI/FPAR within the error range based on a great number of realizations has more statistical significance thus may differ from the specific realization (scene) that was used. Secondly, although the LESS model has been well validated, the confidence of our evaluation results depends on the accuracy of the LESS simulation.
In short, validation in the field of remote sensing utilizing computer simulations has proved feasible. In future studies, we will analyze the other seven biome types, which will provide a more comprehensive evaluation of the MODIS LAI/FPAR retrieval algorithm. In addition, we will change the mode of a single specific scene to obtain retrieval results by simulating multiple scenes. Moreover, evaluation of the algorithm at different levels of vegetation clumping will be the focus of our future research.

Conclusions
This paper presents an uncertainty assessment of the MODIS LAI/FPAR retrieval algorithm over B1 (grassland) based on computer simulation. To accomplish this assessment, we first analyzed the theoretical uncertainty caused by inherent model uncertainty, then we calculated the uncertainty caused by input parameters (BRF and biome type) over simulated 3D grass scenes. Finally, we analyzed the effects of vegetation clumping and scale dependency of the MODIS algorithm. The 3D grass scenes were simulated by a well validated 3D RT model (LESS), which helps to separate the model uncertainty and other uncertainties. We found that the uncertainty of the main and backup algorithm varies considerably. In the same scene, there is a −6.61% bias for the main algorithm retrieval, while the backup algorithm retrieval has a +84.85% bias. We noted that the uncertainty of the saturated retrievals is artificially reduced compared with unsaturated retrievals. At the same time, MODIS showed significant overestimation at low LAI scenes, with a maximum bias of +111.86% for LAI and +162.50% for FPAR. In the high LAI scenes, the "hotspot" geometry results in greater retrieval uncertainty from the backup algorithm. Moreover, input uncertainties further increased the uncertainty of LAI/FPAR retrieval. We found that the uncertainties in BRF in the NIR band has a greater impact than in the red band. The biome type uncertainty also leads to great retrieval uncertainty. Large uncertainties occurred when grassland was misclassified into forest biomes, while smaller uncertainties occurred when the misclassification was within the non-forest biomes. In addition, the clumping effect results in underestimation (−0.846 and −0.525 for the two clumping types, respectively) and we found that the MODIS algorithm is nearly scale-invariant from 100 m to 1000 m pixel sizes. Overall, these results, based on novel computer simulation experiments, can guide the future refinements of the MODIS LAI/FPAR algorithm.