Analysis of Global LAI/FPAR Products from VIIRS and MODIS Sensors for Spatio-Temporal Consistency and Uncertainty from 2012–2016

The operational Moderate Resolution Imaging Spectroradiometer (MODIS) Leaf Area Index (LAI) and Fraction of Photosynthetically Active Radiation absorbed by vegetation (FPAR) algorithm has been successfully implemented for Visible Infrared Imager Radiometer Suite (VIIRS) observations by optimizing a small set of configurable parameters in Look-Up-Tables (LUTs). Our preliminary evaluation showed reasonable agreement between VIIRS and MODIS LAI/FPAR retrievals. However, there is a need for a more comprehensive investigation to assure continuity of multi-sensor global LAI/FPAR time series, as the preliminary evaluation was spatiotemporally limited. In this study, we use a multi-year (2012–2016) global LAI/FPAR product generated from VIIRS and MODIS to evaluate for spatiotemporal consistency. We also quantify uncertainty of the product by utilizing available ground measurements. For both consistency and uncertainty evaluation, we account for variations in biome type and temporal resolution. Our results indicate that the LAI/FPAR retrievals from VIIRS and MODIS are consistent at different spatial (i.e., global and site) and temporal (i.e., 8-day, seasonal and annual) scales. The estimate of mean discrepancy (−0.006 ± 0.013 for LAI and −0.002 ± 0.002 for FPAR) meets the stability requirement for long-term LAI/FPAR Earth System Data Records (ESDRs) from multi-sensors as suggested by the Global Climate Observing System (GCOS). It is noteworthy that the rate of retrievals from the radiative transfer-based main algorithm is also comparable between two sensors. However, a relatively larger discrepancy over tropical forests was observed due to reflectance saturation and an unexpected interannual variation of main algorithm success was noticed due to instability in input surface reflectances. The uncertainties/relative uncertainties of VIIRS and MODIS LAI (FPAR) products assessed through comparisons to ground measurements are estimated to be 0.60/42.2% (0.10/24.4%) and 0.55/39.3% (0.11/26%), respectively. Note that the validated LAI were only distributed in low domains (~2.5), resulting in large relative uncertainty. Therefore, more ground measurements are needed to achieve a more comprehensive evaluation result of product uncertainty. The results presented here generally imbue confidence in the consistency between VIIRS and MODIS LAI/FPAR products and the feasibility of generating long-term multi-sensor LAI/FPAR ESDRs time series.


Introduction
Leaf Area Index (LAI) and Fraction of Photosynthetically Active Radiation absorbed by vegetation (FPAR) are well-known as two essential variables to describe the exchange of fluxes of energy, mass (e.g., water, nutrients, and carbon dioxide) and momentum between the surface and atmosphere [1]. LAI, which is generally defined as one-sided green leaf area per unit ground area in broadleaf canopies and as the projected needle leaf area in coniferous canopies [2], is extensively used to characterize the structure and function of vegetation [3]. FPAR measures the fraction of radiation that leaves absorb in the 0.4-0.7 µm spectrum and it thus evaluates the energy absorption capacity of a canopy [4]. For effective use in most global models of climate, biogeochemistry and ecology, the long-term LAI/FPAR products have been generated from satellite remote sensing [2,5].
In particular, the ground-breaking Earth Observing System (EOS) Moderate Resolution Imaging Spectroradiometer (MODIS) sensors onboard Terra and Aqua satellites have provided an opportunity for opening a new horizon of global LAI/FPAR products [2]. The latest version (Collection 6, C6) of global LAI/FPAR products from MODIS (since February 2000) is freely available through the Land Processes Distributed Active Archive Center (LP DAAC), and is widely used by the scientific, public, and private user communities [5][6][7]. However, both Terra and Aqua MODIS have far exceeded their design life, and they will be likely terminated in the early 2020s [8]. In the context of MODIS termination, the Visible/Infrared Imager Radiometer Suite (VIIRS) instrument onboard the Suomi National Polar-orbiting Partnership (S-NPP) and Joint Polar Satellite System (JPSS) has inherited the scientific roles of MODIS to provide the long-term Earth System Data Records (ESDRs) [9]. Thus, developing consistent LAI/FPAR dataset from these new sensors to continue the MODIS ESDR is a high priority.
Due to sensor-specific spatial resolution and spectral band composition, it is unreliable to directly apply the MODIS algorithm to VIIRS observations [10,11]. With theoretical basis of "canopy spectral invariants" that permits an accurate decoupling of the structural and radiometric components of modeled and measured spectral Bi-directional Reflectance Factors (BRFs) [12][13][14], the MODIS algorithm has been successfully implemented with VIIRS by optimizing a small set of configurable parameters in Look-Up-Tables (LUTs) [11,15]. The preliminary evaluation shows reasonable agreement (global mean difference is 0.024 and 0.029 in January and July of 2015, respectively) between LAI/FPAR retrievals from VIIRS and MODIS [11]. Nevertheless, a comprehensive assessment of the VIIRS product using multi-year and global retrievals is still needed to assure the continuity with the MODIS time series.
Here, we have generated VIIRS LAI/FPAR products that cover the first five-years (2012-2016) of the S-NPP mission era, and we used them to achieve the following objectives: (1) to evaluate the spatiotemporal consistency between the latest version of VIIRS (Version 1, V1) and MODIS (C6) products; (2) to quantify the uncertainty of the VIIRS and MODIS products using available ground measurements; and (3) to understand the observed inconsistency between LAI/FPAR retrievals from two sensors. This paper is organized as follows. Section 2 describes the LAI/FPAR retrieval algorithm and its optimization for VIIRS sensor. Section 3 introduces the data and methods used in this study. Section 4 presents the evaluation results of VIIRS products. Section 5 discusses the reasons for the inconsistency between MODIS and VIIRS products and future developments of VIIRS product evaluation. Finally, Section 6 provides the concluding remarks on this study.

Algorithm Description
The heritage algorithm, developed for MODIS, includes BRFs in the red and Near Infra-Red (NIR) bands, their uncertainties, sun-sensor geometry and a biome classification map. The biome map is required here as a priori-knowledge to reduce the number of unknown parameters of the "ill-posed" inversion problem [2]. The operational LAI/FPAR algorithm includes a main algorithm that is based on a three-dimensional (3D) Radiative Transfer (RT) equation. By describing the photon transfer process, this algorithm links surface spectral BRFs to both structural and spectral parameters of the vegetation canopy and soil [16,17]. Given atmosphere-corrected surface spectral BRFs and their uncertainties, the algorithm finds all candidates of LAI and FPAR by comparing observed and modeled BRFs within biome-specified thresholds of uncertainties that are stored in LUTs. Then, the mean and standard deviation of all LAI/FPAR candidates are reported as the retrieval and its dispersion, respectively. The main algorithm may fail to localize a solution if uncertainties of input BRFs are larger than specified threshold values or due to cloud effects or too low sun/view zenith angles. In this case of a failed solution, the back-up algorithm based on the LAI/FPAR-NDVI relationships for each biome is used to retrieve LAI/FPAR with relatively poor quality [2,13]. Finally, the LAI/FPAR product is generated over the composition period (8-day/4-day) using the daily products. In detail, the compositing algorithm includes two steps: (1) the main algorithm retrievals are first selected (if none available, back-up retrievals are selected); (2) the maximum FPAR and the corresponding LAI are selected as the final product value.

Generation of VIIRS-Specific LUTs
Although VIIRS was designed to be quite similar to MODIS, the differences between these two sensors cannot be ignored [18,19]. For example, the relative spectral responses (RSRs) of MODIS and VIIRS sensors are similar at NIR bands but quite different at red bands (VIIRS has a broader bandwidth than MODIS) [11]. Moreover, VIIRS has a wider swath of 3000 km with a view angle between ±56.28 • than MODIS's swath of 2330 km with a view angle between ±55 • [20], resulting in near-constant resolution of VIIRS from the nadir to the edge of the scan [21]. Additionally, the native instantaneous field of view (IFOV) of VIIRS and MODIS red/NIR bands are 375 m and 250 m, respectively. All discrepancies between VIIRS and MODIS have potential errors for their measured BRFs. Therefore, the generation of consistent LAI/FPAR products from MODIS and VIIRS requires parameterizations that account for sensor-specific features.
The theory of "canopy spectral invariants" provides the necessary parameterizations using a small set of measurable variables that describe the spectral response of vegetation to incident radiation at a range of spatial scales, including: (1) spectrally invariant canopy interceptances, recollision probability, directional escape probability, and their respective hemispherically averaged values; (2) spectrally varying ground reflectance and (3) the single-scattering albedo [12][13][14]. This theory is applied to adapt the LAI/FPAR algorithm to the spatial scale of VIIRS surface reflectance measurements. The spectrally invariant parameters for the different vegetation types in the VIIRS-LUTs were derived using the procedure described in [22]. The spectrally varying ground reflectance was taken from the soil reflectance model developed by [23,24], including 29 typical patterns of effective soil reflectance ranging from bright to dark [15]. Note that the uncertainty of soil reflectance was also considered in generation of VIIRS LUTs to well represent the global conditions [13]. The scale and spectral dependence of the single scattering albedo were exploited to account for VIIRS spatial resolution, spectral bandwidths and radiometric response [25][26][27]. This theory is the fundamental principle of the LAI/FPAR retrieval algorithm because knowing the invariants of the canopy and the single scattering albedo (SSA) of an average phytoelement (e.g., leaf, steam) at any wavelength makes it possible to reconstruct the radiation field of the canopy at any wavelength. The spectral invariant parameters permit decoupling of the structural and radiometric components of any optical sensor signal, which is the theoretical foundation of optimizing configurable parameters to achieve inter-sensor consistency in multi-sensor LAI/FPAR retrievals. Therefore, based on the theory of "canopy spectral invariants", VIIRS-specific LUTs have been developed and incorporated. The theoretical and technical details are described in [11,15].

VIIRS and MODIS LAI/FPAR
Testing and integrating the operational VIIRS LAI/FPAR product are currently under way and its completion is expected in the spring of 2018. The VIIRS LAI/FPAR product will be generated at 500-m spatial resolution and an 8-day time step over a sinusoidal grid. The access to the VIIRS LAI/FPAR product can be obtained from the product's user guide (https://viirsland.gsfc.nasa.gov/PDF/VIIRS_ LAI_FPAR_UserGuide_V1.1.pdf). Here, we used the individual Science Computing Facility (SCF) (Boston University, Boston, MA, USA) to generate a multi-year VIIRS LAI/FPAR product spanning from 2012 to 2016. The operational production system, which is identical to S-NPP Land Science Investigator-led Processing Systems (LandSIPS), was embedded in the SCF, and the optimized LUTs for VIIRS sensor described in Section 2.2 were read for production. The VIIRS Level 2G (L2G) surface reflectance product, called VNP09GA, is composed of all available surface reflectance observations for a given day over a set of tiles with global coverage [28]. By calculating observation scores based on quality assurance (QA) and geometry information (i.e., the basis of high observation coverage, low sensor angle, the absence of cloud, cloud shadow, and aerosol loading), the algorithm produces the intermediate surface reflectance product that contains the best quality observation (i.e., VNP15IP). Then, it accounts for both the best quality observations and the biome map (MCDLCHKM, Figure S1) to produce daily LAI/FPAR retrievals (VNP15A1). With the compositing strategy described in Section 2.1, LAI/FPAR products in the 8-day interval (VNP15A2H) can be obtained. Note that all processes for producing 5-year global VIIRS LAI/FPAR from our SCF are completely identical to the operational production system in S-NPP LandSIPS.
The latest version (C6) of the MODIS LAI/FPAR product was generated and released to the public in 2015. The standard MODIS C6 LAI/FPAR product suite is at 500-m spatial resolution and includes LAI/FPAR retrievals from Terra MODIS, Aqua MODIS, and Terra MODIS + Aqua MODIS Combined [7]. The product suite has two types of temporal compositing periods (i.e., 8 and 4 days) and provides approximately 18-and 16-year records of global land surface from Terra (2000-present) and Aqua (2002-present) MODIS, respectively. The operational production workflow of MODIS is similar to that of the VIIRS product. All the MODIS product suites (MOD15A2H, MYD15A2H, MCD15A2H and MCD15A3H) adopt the same LAI/FPAR retrieval algorithm but different sensors (MOD for Terra and MYD for Aqua) and compositing methods (8-day for A2H and 4-day for A3H). Therefore, LAI/FPAR retrievals for all the MODIS products should be very similar as long as they have comparable observing conditions [29]. To verify this assumption, we compared VIIRS product with different MODIS products and found the results were very similar, that is, different MODIS product suites show similar LAI/FPAR retrievals. As Suomi-NPP is an afternoon-overpassing satellite (13:30 local time), which is identical to MODIS onboard Aqua, we used Aqua MODIS retrieval (MYD15A2H) as a fair counterpart of VIIRS product. Also, the same compositing strategy based on 8-day interval is another reason to use MYD15A2H to compare with the VIIRS LAI/FPAR product. All MODIS LAI/FPAR product suites are available via LP DAAC (https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table) and the details of MYD15A2H can be found in [30].
Both VIIRS and MODIS LAI/FPAR products provide 6 scientific data sets (SDS): (1) Fpar; (2) Lai; (3) FparLai_QC; (4) FparExtra_QC; (5) FparStdDev; and (6) LaiStdDev. The first two SDSs are retrievals of respective FPAR and LAI, and the last two layers are the standard deviation of all candidates from solving the "ill-posed" inversion problem. In addition to LAI/FPAR values, the products contain the corresponding Quality Control (QC) data in the third and fourth layers, and the users are advised to consult the quality flags when using these products [31]. The key indicator of the quality of retrievals is the algorithm path, which distinguishes the following five categories: (1) main algorithm without saturation; (2) main algorithm with saturation; (3) back-up algorithm due to bad geometry (Backup-G); (4) back-up algorithm due to other problems (Backup-O), such as the cloud contamination, snow coverage, etc.; and (5) not produced due to invalid BRFs. In addition to the algorithm path, the QC layers, including the "FparLai_QC" and "FparExtra_QC" data sets, provide information about the presence of clouds, aerosols, and snow, inherited from input reflectance products. Overall, VIIRS and MODIS LAI/FPAR products have very similar characteristics, including the similar input BRFs at red and NIR bands, the similar process chains which consider the good quality of observations by calculating the observation score based on the quality flag (clouds, snow, cloud shadow, etc.) from reflectance products, the same input land cover maps (MCDLCHKM) and the same spatiotemporal resolution (500-m/8-day).

Evaluation of Continuity between VIIRS and MODIS
The evaluation of continuity between VIIRS and MODIS retrievals is prioritized in this study to build a long time series of LAI/FPAR ESDR that is independent of sensor observation. We evaluated the spatiotemporal continuity between VIIRS and MODIS products at two different scales: global and site scale. Continuity is defined as satisfying the following condition, |V−M| < E, where V is the VIIRS estimate of LAI/FPAR, M is the MODIS estimate of LAI/FPAR, and E is an accuracy specification in absolute LAI/FPAR units. In this study, the accuracy specifications for LAI and FPAR are less than 0.25 LAI and 0.02 FPAR, respectively. This specification also meets the long-term ESDR stability requirement suggested by the Global Climate Observing System (GCOS) [32]. However, several issues have to be addressed to achieve a reasonable comparison between VIIRS and MODIS products. First, product geometrical characteristics, including geolocation uncertainties, point spread function (PSF), and projection system, should be considered because the comparison needs to be performed over the same support area [33]. Many studies proposed that the comparison of products over a much larger area than their native spatial resolution is an effective way to reduce potential geolocation errors and PSF impacts [3,34,35]. A good practice suggested by the Committee on Earth Observation Satellites (CEOS) Land Product Validation (LPV) subgroup shows that the product can be evaluated over the area ranged from 3 km × 3 km to 10 km × 10 km [3,[36][37][38]. To understand the product performance over different spatial scales, the VIIRS LAI/FPAR product was evaluated over 10 km × 10 km at global scale and 3 km × 3 km at site scale, respectively. Second, the quality control was performed for VIIRS and MODIS products to exclude pixels contaminated by clouds, shadow, cirrus, and snow [31]. The detailed selection of valid pixels from both VIIRS and MODIS products based on "FparLai_QC" and "FparExtra_QC" layers is shown in Figure S2. Note that only retrievals from the main algorithm were used in our study due to the low accuracy of back-up algorithm retrievals (Section 2.1). For both 10 km × 10 km and 3 km × 3 km areas, the mean LAI/FPAR was calculated only if more than 60% of pixels fell within this area. Additionally, the global land cover map (see Figure S1) was used to identify the vegetated pixel and the specific vegetation type in the analysis of the products.
At the global scale, the aggregated LAI/FPAR retrievals were compared at different temporal scales (8-day, seasonal and annual). We introduced the Benchmark Land Multisite Analysis and Intercomparison of Products (BELMANIP) network, designed to represent the global surface types and conditions, to reduce the geolocation uncertainty and sampling bias [39] at the site scale. The latest version (2.1) of the network contains 445 sites that are located on relatively flat and homogeneous areas within a 10 km × 10 km domain (see Figure S3). Thus, this version is widely used in the product intercomparison practice [31,36]. Similarly, we selected the valid LAI/FPAR retrievals for both VIIRS and MODIS in a 3 km × 3 km area to compare their spatiotemporal consistency using a series of statistical metrics (Bias, Root Mean Square Error (RMSE), relative RMSE (rRMSE) and coefficient of determination (R 2 )). Additionally, the temporal continuity and consistency of the products were also evaluated by biome type. As a key indicator of the quality of retrievals, the algorithm path of each pixel is stored in the LAI/FPAR product. By comparing the retrieval rate of different algorithm paths at the global scale, we can evaluate the overall quality of the product from VIIRS and MODIS. The main algorithm outputs retrievals at high precision in the case of low LAI and at moderate precision when LAI is high and surface reflectance has low sensitivity to LAI. Low-precision retrievals are obtained from the empirical back-up algorithm when the main algorithm fails. An indicator called "Retrieval Index (RI)" that characterizes the proportion of the good quality and high precision retrievals was adopted in our study. The rate of main algorithm retrievals (including retrievals of the main algorithm without saturation, and the main algorithm with saturation) in the 10 km × 10 km area was calculated at the global scale as Equation (1). We investigated RI to show global spatial coverage of the main algorithm retrievals changing with observation dates for the VIIRS product, and also to evaluate whether the proportions of good quality retrievals between VIIRS and MODIS are comparable.

RI =
number o f retrieved pixels via main algorithm number o f total vegetated pixels × 100% (1)

Uncertainty Quantification
The uncertainty is defined as the RMSE between products and ground measurements or ground-measurements-derived reference maps. For this purpose, as most ground measurements were collected during the early EOS era (mostly from 2000 to 2008), here we first synthesized all available in-situ LAI/FPAR measurements over the globe during the VIIRS era (2012-present). Three different field campaign projects were active across 28 sites during the VIIRS era, and those include ImagineS [40], Heihe Watershed Allied Telemetry Experimental (HiWater) [41,42], and Huailai [43]. In particular, the ImagineS project, which is designed to support European Copernicus Global Land Service, has conducted a series of single date or multi-temporal field campaigns over 23 sites around the world during 2013-2016 (http://fp7-imagines.eu). Because of the spatial scale mismatch between ground measurements (tens of meters) and products (500-m), the "bottom-up" approach proposed by the Committee on Earth Observation Satellites (CEOS) Land Product Validation (LPV) subgroup is adopted in this practice [37,44]. The strategy is designed to correlate the scale of in situ biophysical measurements (i.e., LAI and FPAR) to that of the remote sensing product using finer-resolution images to bridge their scale gaps [45]. It is based on a two-stage sampling strategy that (1) uses multiple elementary sampling units (ESUs) to capture the variability across the extent (~3 km × 3 km) of a site and (2) repeats measurements within each ESU to capture the variability within the pixel (~30 m) of high-resolution imagery. Then, an empirical transfer function was established between the ground measurements and the spectral values from the high-resolution imagery to produce LAI/FPAR reference maps. The generated LAI/FPAR maps can be aggregated to match the product pixel grid and are used as the benchmark to validate products [46]. The standard deviation of LAI/FPAR from reference map was also calculated over the 3 km × 3 km region. We also matched the observation date between ground measurements and product at native temporal resolution (8-day). To reduce the geolocation uncertainties caused by different projection systems and point spread functions (PSF), the LAI/FPAR reference maps were averaged over the 3 km ×3 km area for product assessment [33,36]. We collected 68 LAI reference maps and 57 FPAR reference maps from all available ground sites (see Figure S3). However, to reduce unexpected error sources from land surface heterogeneity, we adopted the criterion named information entropy proposed by [31]. The information entropy (E) defined as Equation (2) is to keep only sites with a homogeneous land surface condition.
where ith denotes the specific land cover type from MODIS land cover product, P i indicates the proportion of the area covered by the ith land cover type. The reference map was selected if the E of this site is greater than 1 [31]. Note that the percentage of valid pixels in 3 km × 3 km area should also be greater than 60% as defined in Section 3.2. The detailed information of each reference map is shown in Table S1. Finally, a total of 26 LAI and 24 FPAR reference maps were used in this study. Additionally, the GCOS identified target requirement uncertainties/relative uncertainties of max (0.5, 20%) for LAI and max (0.05, 10%) for FPAR [32], which will also be utilized as the criterion to evaluate the uncertainty of VIIRS and MODIS products.

Spatiotemporal Consistency between VIIRS and MODIS
4.1.1. Global Scale Figure 1 depicts the 8-day global LAI/FPAR time series of VIIRS and MODIS from 2012 to 2016. Both VIIRS and MODIS tightly well capture the LAI and FPAR seasonality, with two distinct peaks every year caused by the hemispherically different seasons. The LAI and FPAR in boreal summer (i.e., July) are much greater (about 2.0 LAI and 0.51 FPAR) than those (about 1.5 LAI and 0.42 FPAR) in boreal winter (i.e., January) because of larger proportion of global vegetation in northern hemisphere. In particular, both retrievals from VIIRS and MODIS agree very well with long-term time series, with a mean difference of −0.006 ± 0.013 for LAI and −0.002 ± 0.002 for FPAR. The LAI/FPAR difference between VIIRS and MODIS meets the stability requirement of multi-sensor product time-series suggested by GCOS [32]. The standard deviation (std.) of the difference illustrates a slight stochastic difference between VIIRS and MODIS products. The seasonal change of std. of difference is also observed with higher values (0.46 for LAI and 0.060 for FPAR) in January and lower values (0.36 for LAI and 0.045 for FPAR) in July. This is because in January, the LAI/FPAR differences are lower in most regions but larger in tropical areas, resulting in larger std. of difference, as shown in Figure 2. (0.5, 20%) for LAI and max (0.05, 10%) for FPAR [32], which will also be utilized as the criterion to evaluate the uncertainty of VIIRS and MODIS products.

Spatiotemporal Consistency between VIIRS and MODIS
4.1.1. Global Scale Figure 1 depicts the 8-day global LAI/FPAR time series of VIIRS and MODIS from 2012 to 2016. Both VIIRS and MODIS tightly well capture the LAI and FPAR seasonality, with two distinct peaks every year caused by the hemispherically different seasons. The LAI and FPAR in boreal summer (i.e., July) are much greater (about 2.0 LAI and 0.51 FPAR) than those (about 1.5 LAI and 0.42 FPAR) in boreal winter (i.e., January) because of larger proportion of global vegetation in northern hemisphere. In particular, both retrievals from VIIRS and MODIS agree very well with long-term time series, with a mean difference of −0.006 ± 0.013 for LAI and −0.002 ± 0.002 for FPAR. The LAI/FPAR difference between VIIRS and MODIS meets the stability requirement of multi-sensor product time-series suggested by GCOS [32]. The standard deviation (std.) of the difference illustrates a slight stochastic difference between VIIRS and MODIS products. The seasonal change of std. of difference is also observed with higher values (0.46 for LAI and 0.060 for FPAR) in January and lower values (0.36 for LAI and 0.045 for FPAR) in July. This is because in January, the LAI/FPAR differences are lower in most regions but larger in tropical areas, resulting in larger std. of difference, as shown in Figure 2. The top panels of (a,b) show the seasonal variation of VIIRS and MODIS LAI/FPAR retrievals, while the bottom two panels show the mean and standard deviation (std.) of difference (the difference was first calculated in each 10 km × 10 km grid, and then the mean and std. were calculated for all grids) between the two products (VIIRS-MODIS). The mean difference indicates whether VIIRS and MODIS products have systematic errors, while the std. of difference shows the dispersion of errors between VIIRS and MODIS products. Also see Figure S4 for the evaluation between VIIRS and MODIS LAI/FPAR products at the native spatial resolution (500-m) for each biome.
Details of the comparison at the seasonal scale can be found in Table 1. The mean and standard deviation (std.) of product difference were calculated at a 3-month interval (March-May (MAM), June-August (JJA), September-November (SON) and December-February (DJF)). Both LAI and FPAR show the largest differences in JJA, which are primarily due to the largest LAI and FPAR values in this period ( Figure 1). Nevertheless, the std. presents the lowest magnitude in JJA, showing the similar differences at the global scale. Although the differences between VIIRS and MODIS change in different seasons, the differences present low mean and std. values. Results show that the LAI/FPAR difference between products can also satisfy the stability requirement of a The top panels of (a,b) show the seasonal variation of VIIRS and MODIS LAI/FPAR retrievals, while the bottom two panels show the mean and standard deviation (std.) of difference (the difference was first calculated in each 10 km × 10 km grid, and then the mean and std. were calculated for all grids) between the two products (VIIRS-MODIS). The mean difference indicates whether VIIRS and MODIS products have systematic errors, while the std. of difference shows the dispersion of errors between VIIRS and MODIS products. Also see Figure S4 for the evaluation between VIIRS and MODIS LAI/FPAR products at the native spatial resolution (500-m) for each biome.
Details of the comparison at the seasonal scale can be found in Table 1. The mean and standard deviation (std.) of product difference were calculated at a 3-month interval (March-May (MAM), June-August (JJA), September-November (SON) and December-February (DJF)). Both LAI and FPAR show the largest differences in JJA, which are primarily due to the largest LAI and FPAR values in this  (Figure 1). Nevertheless, the std. presents the lowest magnitude in JJA, showing the similar differences at the global scale. Although the differences between VIIRS and MODIS change in different seasons, the differences present low mean and std. values. Results show that the LAI/FPAR difference between products can also satisfy the stability requirement of a long-term dataset for different seasons. Additionally, the mean and std. of LAI/FPAR difference also depend on the biome type (see Table S2). From biome-specific investigations, we do not find obvious systematic over-or under-estimation across biomes. Forest biomes (B5-B8) present relatively larger LAI/FPAR differences than non-forest biomes (B1-B4) in each season. Among all biomes, the needleleaf forest (B7-B8) shows the largest absolute mean difference between products, while the evergreen broadleaf forest (B5) presents the largest std. of LAI/FPAR difference. As shown in Table 1, most biomes show the largest mean difference in JJA among four seasons except for the needleleaf forest (B7-B8), which shows the largest differences in DJF. This is likely due to the snow coverage and the large solar zenith angle in the high latitude areas in DJF, introducing larger uncertainty in surface reflectances. The possible reasons for LAI/FPAR inconsistencies between VIIRS and MODIS products will be further discussed in Section 5.1. The "MAM", "JJA", "SON" and "DJF" stand for "March-May", "June-August", "September-November" and "December-February", respectively.
In detail, the global spatial patterns of mean LAI and FPAR difference maps between VIIRS and MODIS products at two typical observation dates (17-24 January and 12-19 July) for the five-year period (2012-2016) are presented in Figure 2. The detailed spatial distribution of LAI/FPAR differences shows no observed systematic bias between VIIRS and MODIS products. The overall differences of LAI and FPAR are within ±0.25 (percentage: 83%) and ±0.02 (percentage: 70%), respectively, with larger differences over the dense forests. The invalid LAI/FPAR values at high latitudes in the northern hemisphere in January are mainly due to the snow coverage, polar night, and large solar zenith angle [11]. For LAI retrievals (Figure 2a,b), most regions in the northern hemisphere have a slightly larger difference in July than in January. However, the LAI difference in the equatorial region, e.g., Amazonia and central Africa, shows an opposite pattern due to the growing season of vegetation in January. From the spatial pattern comparison of LAI differences in January and July at the global scale, LAI differences depend on the LAI magnitudes, with larger LAI differences in high LAI domains. The possible reasons will be discussed in Section 5.1. We also observed similar spatial distribution of FPAR differences from Figure 2c,d. long-term dataset for different seasons. Additionally, the mean and std. of LAI/FPAR difference also depend on the biome type (see Table S2). From biome-specific investigations, we do not find obvious systematic over-or under-estimation across biomes. Forest biomes (B5-B8) present relatively larger LAI/FPAR differences than non-forest biomes (B1-B4) in each season. Among all biomes, the needleleaf forest (B7-B8) shows the largest absolute mean difference between products, while the evergreen broadleaf forest (B5) presents the largest std. of LAI/FPAR difference. As shown in Table 1, most biomes show the largest mean difference in JJA among four seasons except for the needleleaf forest (B7-B8), which shows the largest differences in DJF. This is likely due to the snow coverage and the large solar zenith angle in the high latitude areas in DJF, introducing larger uncertainty in surface reflectances. The possible reasons for LAI/FPAR inconsistencies between VIIRS and MODIS products will be further discussed in Section 5.1. The "MAM", "JJA", "SON" and "DJF" stand for "March-May", "June-August", "September-November" and "December-February", respectively.
In detail, the global spatial patterns of mean LAI and FPAR difference maps between VIIRS and MODIS products at two typical observation dates (17-24 January and 12-19 July) for the five-year period (2012-2016) are presented in Figure 2. The detailed spatial distribution of LAI/FPAR differences shows no observed systematic bias between VIIRS and MODIS products. The overall differences of LAI and FPAR are within ±0.25 (percentage: 83%) and ±0.02 (percentage: 70%), respectively, with larger differences over the dense forests. The invalid LAI/FPAR values at high latitudes in the northern hemisphere in January are mainly due to the snow coverage, polar night, and large solar zenith angle [11]. For LAI retrievals (Figure 2a,b), most regions in the northern hemisphere have a slightly larger difference in July than in January. However, the LAI difference in the equatorial region, e.g., Amazonia and central Africa, shows an opposite pattern due to the growing season of vegetation in January. From the spatial pattern comparison of LAI differences in January and July at the global scale, LAI differences depend on the LAI magnitudes, with larger LAI differences in high LAI domains. The possible reasons will be discussed in Section 5.1. We also observed similar spatial distribution of FPAR differences from Figure 2c,d. As VIIRS will continue the role of MODIS to provide LAI/FPAR records, it is necessary to evaluate whether VIIRS can inherit MODIS for the long-term analysis of vegetation dynamics (i.e., LAI trend over several years [47,48]) at global scale. Here, we additionally investigate how well VIIRS captures the interannual variability of MODIS retrievals by splitting three large regions: the Northern Hemisphere (25° N-90° N), the equatorial region (25° N-25° S) and the Southern Hemisphere (25° S-90° S). The LAI anomaly, defined as a certain year of mean LAI value minus mean LAI value over 2012-2016, is used to quantify the interannual variability of products. A large LAI anomaly indicates that LAI changes considerably compared with the mean LAI over different years. As shown in Figure 3, VIIRS shows quite similar interannual variation for each region to those of MODIS. However, the magnitude of LAI anomalies is slightly different, especially in the equatorial region (Figure 3b). The variation (~0.03) of MODIS LAI anomaly in the equatorial region is much larger than that for VIIRS (~0.01) from 2013 to 2014, resulting in the opposite variation of 2014 at global scale (Figure 3d). This is likely because VIIRS reflectance is easier to be saturated due to the sensor's intrinsic limitation (VIIRS reflectance is generally larger in the NIR band and lower in the red band [11]). The saturated reflectance over dense canopy in tropical forests can only provide limited information for LAI retrievals [49], that is, the real vegetation change cannot be shown if the reflectance is saturated, resulting in the small variation of VIIRS LAI retrievals over different years. As VIIRS will continue the role of MODIS to provide LAI/FPAR records, it is necessary to evaluate whether VIIRS can inherit MODIS for the long-term analysis of vegetation dynamics (i.e., LAI trend over several years [47,48]) at global scale. Here, we additionally investigate how well VIIRS captures the interannual variability of MODIS retrievals by splitting three large regions: the Northern Hemisphere (25 • N-90 • N), the equatorial region (25 • N-25 • S) and the Southern Hemisphere (25 • S-90 • S). The LAI anomaly, defined as a certain year of mean LAI value minus mean LAI value over 2012-2016, is used to quantify the interannual variability of products. A large LAI anomaly indicates that LAI changes considerably compared with the mean LAI over different years. As shown in Figure 3, VIIRS shows quite similar interannual variation for each region to those of MODIS. However, the magnitude of LAI anomalies is slightly different, especially in the equatorial region (Figure 3b). The variation (~0.03) of MODIS LAI anomaly in the equatorial region is much larger than that for VIIRS (~0.01) from 2013 to 2014, resulting in the opposite variation of 2014 at global scale (Figure 3d). This is likely because VIIRS reflectance is easier to be saturated due to the sensor's intrinsic limitation (VIIRS reflectance is generally larger in the NIR band and lower in the red band [11]). The saturated reflectance over dense canopy in tropical forests can only provide limited information for LAI retrievals [49], that is, the real vegetation change cannot be shown if the reflectance is saturated, resulting in the small variation of VIIRS LAI retrievals over different years. As VIIRS will continue the role of MODIS to provide LAI/FPAR records, it is necessary to evaluate whether VIIRS can inherit MODIS for the long-term analysis of vegetation dynamics (i.e., LAI trend over several years [47,48]) at global scale. Here, we additionally investigate how well VIIRS captures the interannual variability of MODIS retrievals by splitting three large regions: the Northern Hemisphere (25° N-90° N), the equatorial region (25° N-25° S) and the Southern Hemisphere (25° S-90° S). The LAI anomaly, defined as a certain year of mean LAI value minus mean LAI value over 2012-2016, is used to quantify the interannual variability of products. A large LAI anomaly indicates that LAI changes considerably compared with the mean LAI over different years. As shown in Figure 3, VIIRS shows quite similar interannual variation for each region to those of MODIS. However, the magnitude of LAI anomalies is slightly different, especially in the equatorial region (Figure 3b). The variation (~0.03) of MODIS LAI anomaly in the equatorial region is much larger than that for VIIRS (~0.01) from 2013 to 2014, resulting in the opposite variation of 2014 at global scale (Figure 3d). This is likely because VIIRS reflectance is easier to be saturated due to the sensor's intrinsic limitation (VIIRS reflectance is generally larger in the NIR band and lower in the red band [11]). The saturated reflectance over dense canopy in tropical forests can only provide limited information for LAI retrievals [49], that is, the real vegetation change cannot be shown if the reflectance is saturated, resulting in the small variation of VIIRS LAI retrievals over different years.   (Figure 4), non-forest biomes (i.e., B1-B4) have better agreement between products than forest biomes, which can also be found for FPAR retrievals ( Figure 5). The best agreement is achieved over shrubs (Bias = 0.009, RMSE = 0.060) and broadleaf crops (Bias = −0.001, RMSE = 0.029) for LAI and FPAR retrievals, respectively. VIIRS retrievals from most biomes, except evergreen broadleaf forest (B5), capture more than 85% of the variations (R 2 > 0.85) observed in MODIS retrievals. The case of dense evergreen broadleaf forest shows less agreement in both LAI (R 2 = 0.62, RMSE = 0.69) and FPAR (R 2 = 0.47, RMSE = 0.048), which is not surprising because of sub-optimal quality of retrievals due to reflectance saturation for dense forest.
In detail, small variations in reflectances can result in large variation in LAI under the condition of saturation [49]. Additionally, another possible reason of the inconsistency between VIIRS and MODIS retrievals for evergreen broadleaf forest will be discussed in Section 5.1.

Site Scale
The VIIRS and MODIS products are further compared at the site scale for each biome over all BELMANIP-2.1 sampling sites. Figures 4 and 5 show the density scatter plots between VIIRS and MODIS for all 8-day LAI and FPAR retrievals from 2012 to 2016, respectively. The statistical results (R 2 , Bias and RMSE) of all pixels indicate good consistency between VIIRS and MODIS in both LAI/FPAR retrievals. For LAI retrievals (Figure 4), non-forest biomes (i.e., B1-B4) have better agreement between products than forest biomes, which can also be found for FPAR retrievals ( Figure 5). The best agreement is achieved over shrubs (Bias = 0.009, RMSE = 0.060) and broadleaf crops (Bias = −0.001, RMSE = 0.029) for LAI and FPAR retrievals, respectively. VIIRS retrievals from most biomes, except evergreen broadleaf forest (B5), capture more than 85% of the variations (R 2 > 0.85) observed in MODIS retrievals. The case of dense evergreen broadleaf forest shows less agreement in both LAI (R 2 = 0.62, RMSE = 0.69) and FPAR (R 2 = 0.47, RMSE = 0.048), which is not surprising because of sub-optimal quality of retrievals due to reflectance saturation for dense forest. In detail, small variations in reflectances can result in large variation in LAI under the condition of saturation [49]. Additionally, another possible reason of the inconsistency between VIIRS and MODIS retrievals for evergreen broadleaf forest will be discussed in Section 5.1.

Site Scale
The VIIRS and MODIS products are further compared at the site scale for each biome over all BELMANIP-2.1 sampling sites. Figures 4 and 5 show the density scatter plots between VIIRS and MODIS for all 8-day LAI and FPAR retrievals from 2012 to 2016, respectively. The statistical results (R 2 , Bias and RMSE) of all pixels indicate good consistency between VIIRS and MODIS in both LAI/FPAR retrievals. For LAI retrievals (Figure 4), non-forest biomes (i.e., B1-B4) have better agreement between products than forest biomes, which can also be found for FPAR retrievals ( Figure 5). The best agreement is achieved over shrubs (Bias = 0.009, RMSE = 0.060) and broadleaf crops (Bias = −0.001, RMSE = 0.029) for LAI and FPAR retrievals, respectively. VIIRS retrievals from most biomes, except evergreen broadleaf forest (B5), capture more than 85% of the variations (R 2 > 0.85) observed in MODIS retrievals. The case of dense evergreen broadleaf forest shows less agreement in both LAI (R 2 = 0.62, RMSE = 0.69) and FPAR (R 2 = 0.47, RMSE = 0.048), which is not surprising because of sub-optimal quality of retrievals due to reflectance saturation for dense forest.
In detail, small variations in reflectances can result in large variation in LAI under the condition of saturation [49]. Additionally, another possible reason of the inconsistency between VIIRS and MODIS retrievals for evergreen broadleaf forest will be discussed in Section 5.1.        Figure 6 shows the seasonal trajectory of VIIRS and MODIS LAI/FPAR over four example sites (grasses/cereal crops (B1), shrubs (B2), deciduous broadleaf forest (B6) and deciduous needleleaf forest (B8)) for the period from 2012 to 2016. All available ground measurements are additionally plotted as a reference. Both products shown in Figure 6 achieve a good temporal consistency across all biomes and comparable spatial coverage (i.e., missing data rate). Each site shows distinct characteristics of seasonal variations. For example, non-forest biomes (Figure 6a,b) exhibit smaller seasonal amplitude, less than 2 for LAI and 0.4 for FPAR, while forest biomes (Figure 6c,d) reveal a higher LAI (about 4 LAI unit) and FPAR (about 0.5 FPAR unit) seasonality. In addition, the algorithm implemented for both sensors agrees quite well with the ground measurements, indicating VIIRS and MODIS can capture the temporal trajectory of true LAI/FPAR values.   Figure 6 shows the seasonal trajectory of VIIRS and MODIS LAI/FPAR over four example sites (grasses/cereal crops (B1), shrubs (B2), deciduous broadleaf forest (B6) and deciduous needleleaf forest (B8)) for the period from 2012 to 2016. All available ground measurements are additionally plotted as a reference. Both products shown in Figure 6 achieve a good temporal consistency across all biomes and comparable spatial coverage (i.e., missing data rate). Each site shows distinct characteristics of seasonal variations. For example, non-forest biomes (Figure 6a,b) exhibit smaller seasonal amplitude, less than 2 for LAI and 0.4 for FPAR, while forest biomes (Figure 6c,d) reveal a higher LAI (about 4 LAI unit) and FPAR (about 0.5 FPAR unit) seasonality. In addition, the algorithm implemented for both sensors agrees quite well with the ground measurements, indicating VIIRS and MODIS can capture the temporal trajectory of true LAI/FPAR values.  Figure 6 shows the seasonal trajectory of VIIRS and MODIS LAI/FPAR over four example sites (grasses/cereal crops (B1), shrubs (B2), deciduous broadleaf forest (B6) and deciduous needleleaf forest (B8)) for the period from 2012 to 2016. All available ground measurements are additionally plotted as a reference. Both products shown in Figure 6 achieve a good temporal consistency across all biomes and comparable spatial coverage (i.e., missing data rate). Each site shows distinct characteristics of seasonal variations. For example, non-forest biomes (Figure 6a,b) exhibit smaller seasonal amplitude, less than 2 for LAI and 0.4 for FPAR, while forest biomes (Figure 6c,d) reveal a higher LAI (about 4 LAI unit) and FPAR (about 0.5 FPAR unit) seasonality. In addition, the algorithm implemented for both sensors agrees quite well with the ground measurements, indicating VIIRS and MODIS can capture the temporal trajectory of true LAI/FPAR values.

Spatial Coverage
The VIIRS RI and the Difference of RI (DRI) between VIIRS and MODIS were calculated at the global scale from 2012 to 2016 to ensure their consistency of spatial coverage. For the sake of brevity, here we only present the results from two 8-day compositional dates (i.e., 17-24 January and 12-19 July) representing the respective boreal winter and summer season. Figure 7 shows seasonally different biome-specific retrieval success (i.e., RI) and RI difference (DRI) between VIIRS and MODIS. The global RIs of VIIRS in January and July are 58.8% and 85.1%, respectively. In general, the RI of forest biomes (B5-B8) is lower than that of non-forest biomes (B1-B4) for both two seasons. Both sensors give a lower RI in January for most biomes (see Figure S5), especially for needleleaf forest (B7-B8, <20%). The lower RI is partly because most needleleaf forests are situated in the high latitude or altitude regions (see Figure S1) where they are often covered by snow in January. In addition, a larger solar zenith angle at high latitudes in January is the most important reason for low needleleaf forest RI (See Figure S5c,d). The evergreen broadleaf forests (B5) in the equatorial region have comparable RIs in January and July, and relatively lower RIs than other biomes due to the impact of cloudy condition throughout the year.
From the comparison with calculated MODIS RIs, we find that VIIRS completely resembles their seasonal variations of each algorithm pathway including main and backup (Backup-G and Backup-O) retrievals (see details in Figure S5). As shown in Figure 7, the proportion of DRI ranging from −5% to 5% at global scale is 84.7% in January and 77.2% in July. For individual biome types, both VIIRS and MODIS RIs agree well in January for each biome type showing a low mean DRI (~1%). However, we also observed a relatively large mean DRI (~5%) from the comparison in July, especially for forest biomes. The reason for this inconsistency of RIs will be discussed in Section 5.1. The standard deviation (std.) of DRI shows that VIIRS has a stable DRI (~4%) with MODIS for shrubs (B2), while also presents a relatively unstable DRI (~9%) for evergreen broadleaf forest (B5) over five years.

Spatial Coverage
The VIIRS RI and the Difference of RI (DRI) between VIIRS and MODIS were calculated at the global scale from 2012 to 2016 to ensure their consistency of spatial coverage. For the sake of brevity, here we only present the results from two 8-day compositional dates (i.e., 17-24 January and 12-19 July) representing the respective boreal winter and summer season. Figure 7 shows seasonally different biome-specific retrieval success (i.e., RI) and RI difference (DRI) between VIIRS and MODIS. The global RIs of VIIRS in January and July are 58.8% and 85.1%, respectively. In general, the RI of forest biomes (B5-B8) is lower than that of non-forest biomes (B1-B4) for both two seasons. Both sensors give a lower RI in January for most biomes (see Figure S5), especially for needleleaf forest (B7-B8, <20%). The lower RI is partly because most needleleaf forests are situated in the high latitude or altitude regions (see Figure S1) where they are often covered by snow in January. In addition, a larger solar zenith angle at high latitudes in January is the most important reason for low needleleaf forest RI (See Figure S5c,d). The evergreen broadleaf forests (B5) in the equatorial region have comparable RIs in January and July, and relatively lower RIs than other biomes due to the impact of cloudy condition throughout the year.
From the comparison with calculated MODIS RIs, we find that VIIRS completely resembles their seasonal variations of each algorithm pathway including main and backup (Backup-G and Backup-O) retrievals (see details in Figure S5). As shown in Figure 7, the proportion of DRI ranging from −5% to 5% at global scale is 84.7% in January and 77.2% in July. For individual biome types, both VIIRS and MODIS RIs agree well in January for each biome type showing a low mean DRI (~1%). However, we also observed a relatively large mean DRI (~5%) from the comparison in July, especially for forest biomes. The reason for this inconsistency of RIs will be discussed in Section 5.1. The standard deviation (std.) of DRI shows that VIIRS has a stable DRI (~4%) with MODIS for shrubs (B2), while also presents a relatively unstable DRI (~9%) for evergreen broadleaf forest (B5) over five years.

Uncertainty Assessment
The uncertainty of both VIIRS and MODIS LAI/FPAR products was evaluated using the synthesized ground measurements, as shown in Figure 8. In general, our comparison reveals that VIIRS and MODIS show comparable uncertainties of LAI (VIIRS = 0.60, MODIS = 0.55) and FPAR (VIIRS = 0.10, MODIS = 0.11) retrievals. As the result of MODIS resembles that of VIIRS, here, we use

Uncertainty Assessment
The uncertainty of both VIIRS and MODIS LAI/FPAR products was evaluated using the synthesized ground measurements, as shown in Figure 8. In general, our comparison reveals that VIIRS and MODIS show comparable uncertainties of LAI (VIIRS = 0.60, MODIS = 0.55) and FPAR (VIIRS = 0.10, MODIS = 0.11) retrievals. As the result of MODIS resembles that of VIIRS, here, we use VIIRS as a case to explain the general pattern of the results. For LAI, both ground measurements corrected (hereafter, true LAI) or uncorrected (hereafter, effective LAI) for clumping effect are introduced in this study but separately analyzed. Since the clumping effect was taken into account in the LAI retrieval algorithm, VIIRS displays better agreement with true LAI (Bias = 0.17, RMSE = 0.60, rRMSE = 42.2%) than with effective LAI (bias = 0.47, RMSE = 0.79, rRMSE = 66.0%) as expected. Specifically, VIIRS shows an important overestimation tendency for the effective LAI case (Bias = 0.47) because the ground-measured effective LAI is generally lower than true LAI [50]. As LAI was only distributed on relatively lower domains (<2.5), the RMSE is better than the rRMSE to show LAI product accuracy based on the accuracy requirement (max(0.5, 20%)) from GCOS. For FPAR, most scatters are distributed within 0-0.2 difference, indicating an overestimation tendency of VIIRS FPAR (Bias = 0.08). Nevertheless, VIIRS FPAR retrieval displays good agreement with ground measured FPAR in terms of R 2 (=0.83), RMSE (=0.10) and rRMSE (=24.4%), indicating that VIIRS can well capture the FPAR variation (the dispersion of VIIRS FPAR is low). In addition, we can see that both VIIRS and MODIS products cannot totally meet the target accuracy requirements suggested by GCOS, with 66.7% LAI and 37.5% FPAR pixels less than the LAI (max (0.5, 20%)) and FPAR (max (0.05, 10%)) specific criteria, respectively. However, it should be noted that the uncertainty of LAI/FPAR retrievals comes from both VIIRS products and other sources including uncertainties of reference maps and a mismatch in the spatial and temporal domains, despite that efforts were made to reduce such sources of uncertainty. Additionally, we should also note that the distribution of measurements is not sufficient enough because only small number of measurements and no forest biomes are included. Adding more ground measurements, especially from forest biomes, may achieve a more robust evaluation result of VIIRS LAI/FPAR product. Further discussions will be followed in Section 5.2. VIIRS as a case to explain the general pattern of the results. For LAI, both ground measurements corrected (hereafter, true LAI) or uncorrected (hereafter, effective LAI) for clumping effect are introduced in this study but separately analyzed. Since the clumping effect was taken into account in the LAI retrieval algorithm, VIIRS displays better agreement with true LAI (Bias = 0.17, RMSE = 0.60, rRMSE = 42.2%) than with effective LAI (bias = 0.47, RMSE = 0.79, rRMSE = 66.0%) as expected. Specifically, VIIRS shows an important overestimation tendency for the effective LAI case (Bias = 0.47) because the ground-measured effective LAI is generally lower than true LAI [50]. As LAI was only distributed on relatively lower domains (<2.5), the RMSE is better than the rRMSE to show LAI product accuracy based on the accuracy requirement (max(0.5, 20%)) from GCOS. For FPAR, most scatters are distributed within 0-0.2 difference, indicating an overestimation tendency of VIIRS FPAR (Bias = 0.08). Nevertheless, VIIRS FPAR retrieval displays good agreement with ground measured FPAR in terms of R 2 (=0.83), RMSE (=0.10) and rRMSE (=24.4%), indicating that VIIRS can well capture the FPAR variation (the dispersion of VIIRS FPAR is low). In addition, we can see that both VIIRS and MODIS products cannot totally meet the target accuracy requirements suggested by GCOS, with 66.7% LAI and 37.5% FPAR pixels less than the LAI (max (0.5, 20%)) and FPAR (max (0.05, 10%)) specific criteria, respectively. However, it should be noted that the uncertainty of LAI/FPAR retrievals comes from both VIIRS products and other sources including uncertainties of reference maps and a mismatch in the spatial and temporal domains, despite that efforts were made to reduce such sources of uncertainty. Additionally, we should also note that the distribution of measurements is not sufficient enough because only small number of measurements and no forest biomes are included. Adding more ground measurements, especially from forest biomes, may achieve a more robust evaluation result of VIIRS LAI/FPAR product. Further discussions will be followed in Section 5.2.  The red dashed lines show the GCOS specifications boundaries (max (0.5, 20%) for LAI and max (0.05, 10%) for FPAR). Note that LAI was only distributed in relatively lower domains (~2.5), the RMSE is better than the rRMSE to show LAI product accuracy based on the accuracy requirement (max (0.5, 20%)) from GCOS.

Understanding Inconsistency between VIIRS and MODIS
From our comprehensive evaluation presented in Section 4, we can conclude that the implemented LAI/FPAR algorithm on the VIIRS can generate highly consistent LAI/FPAR retrievals with those of MODIS at multi-spatiotemporal scales. However, as we observed some minor discrepancies between VIIRS and MODIS products; here, we investigate the possible reasons to better understand observed inconsistencies. The saturation frequency and dispersion of the retrieved LAI distribution (DLAI) are two elements by which the quality of the retrieval can be assessed. The accuracy of the retrieval decreases under conditions of saturation, that is, the reflectance data do not contain accurate information about the surface [13,49]. DLAI is defined as the standard deviation of all LAI candidates in LUTs enabling to reproduce input surface reflectances under various canopy and soil combinations (recall that the algorithm is solving the ill-posed inverse problem) [13,15]. Thus, in theory, the lower the DLAI value, the higher the accuracy of the algorithm as DLAI represents the degree of difficulty to localize a single solution.
In this manner, Figure 9 shows important relationships between observed difference and reliability of retrievals (i.e., saturation rate and DLAI). Note that all quantities in Figure 9 are based on 8-day global products from 2012 to 2016. First, Figure 9a shows how the saturation rate behaves as a function of LAI. Unsurprisingly, the saturation rate increases with increasing LAI and thus we can infer the accuracy of the retrievals is decreasing. In the same vein, DLAI increases with increasing LAI but DLAI tend to decrease at certain LAI (about 5 for non-forest and 4.5 for forest). This is because DLAI is not only related to LAI (higher LAI yields higher DLAI), but also depends on the number of available solutions in LUTs (Section 2.1), i.e., the nature of the standard deviation. Under the relationship of saturation rate and DLAI to LAI retrieval, we found that observed discrepancy between VIIRS and MODIS tends to increase with an increasing saturation rate and DLAI. We highlight that the discrepancy is abruptly increased when the saturation rate shows a sharp rise in higher LAI, namely, the reliability of retrieval determines the degree of consistency between two sensors. Additionally, we should also note the LAI difference (Bias = 0.02, RMSE = 0.69) for evergreen broadleaf forest (EBF) is not negligible (see Figure 4e); thus, a special caution is needed for this biome. Observing land surface over tropical regions where EBF mostly situated is a classical challenge in optical remote sensing field due to cloudy or high aerosol atmospheric conditions [51,52]. Therefore, this difference we observed here can be likely explained by different observing conditions in two sensors. For instance, although both Suomi-NPP and Aqua are afternoon-overpassing satellites, different orbital path and selection of daily best input observation may result in a dissimilar input surface reflectance. The impact of such discrepancy is expected to increase in highly varying atmospheric conditions over densely vegetated tropical forests. Another possible reason is QA which is slightly different between VIIRS and MODIS. For instance, Platnick et al. (2016) reported that cloud detection and flagging its mask in VIIRS is limited by absence of CO 2 and H 2 O absorption channels [53], in turn such discrepancy results in unlike LAI/FPAR retrievals.
In addition, we also observed the relatively large RI difference in July (see Figure 7), which is likely explained by the reported instability of input surface reflectance for VIIRS [8]. The results from [8] show larger deviations of spectral adjustment coefficients to make MODIS-like VIIRS surface reflectances are observed for the initial years of VIIRS operation, while relatively better stability is observed for later years. As VIIRS LUTs have been developed using BRFs in 2015 based on the time-invariant basis, unstable input surface reflectances over different years may induce unexpected interannual variability of RIs. Therefore, the DRI between VIIRS and MODIS is less than 1% in 2015, whereas it shows relatively larger difference during 2012-2014 (see Figure S6). Here, we should note that RI differences of January and non-forest are smaller than July and forest biomes ( Figure S6). The lower RI values are caused by the algorithm path, which is also sensitive to the high LAI values. More specifically, the domain of the main algorithm retrievals is more limited in dense vegetation with higher NIR reflectance and lower red reflectance (see Figure 10a in [7]), which introduces lower RI due to the interannual variability of surface reflectance.
There are additional minor factors that influence the inter-consistency of the LAI/FPAR product: compositing and off-nadir view zenith angle (VZA), which are from the intrinsic limitation of VIIRS and MODIS sensors. For both VIIRS and MODIS, the algorithm first accounts for 8 daily surface reflectances and produces daily LAI/FPAR. Based on the maximum FPAR compositing strategy, a daily retrieval is selected as 8-day product value. When selected VIIRS daily retrieval is separate from MODIS daily retrieval, the degree of inconsistency between VIIRS and MODIS may increase as [36] shows similar pattern in surface reflectance. For the off-nadir VZA case, MODIS pixels increase by a factor of six from nadir to edge of scan, whereas VIIRS restricts the pixel to an approximate twofold increase using an onboard aggregation scheme [11]. Therefore, a reduced effective spatial resolution of MODIS with the increase of VZA may introduce additional source of inconsistency between VIIRS and MODIS products [54].
on the time-invariant basis, unstable input surface reflectances over different years may induce unexpected interannual variability of RIs. Therefore, the DRI between VIIRS and MODIS is less than 1% in 2015, whereas it shows relatively larger difference during 2012-2014 (see Figure S6). Here, we should note that RI differences of January and non-forest are smaller than July and forest biomes ( Figure S6). The lower RI values are caused by the algorithm path, which is also sensitive to the high LAI values. More specifically, the domain of the main algorithm retrievals is more limited in dense vegetation with higher NIR reflectance and lower red reflectance (see Figure 10a in [7]), which introduces lower RI due to the interannual variability of surface reflectance.
There are additional minor factors that influence the inter-consistency of the LAI/FPAR product: compositing and off-nadir view zenith angle (VZA), which are from the intrinsic limitation of VIIRS and MODIS sensors. For both VIIRS and MODIS, the algorithm first accounts for 8 daily surface reflectances and produces daily LAI/FPAR. Based on the maximum FPAR compositing strategy, a daily retrieval is selected as 8-day product value. When selected VIIRS daily retrieval is separate from MODIS daily retrieval, the degree of inconsistency between VIIRS and MODIS may increase as [36] shows similar pattern in surface reflectance. For the off-nadir VZA case, MODIS pixels increase by a factor of six from nadir to edge of scan, whereas VIIRS restricts the pixel to an approximate twofold increase using an onboard aggregation scheme [11]. Therefore, a reduced effective spatial resolution of MODIS with the increase of VZA may introduce additional source of inconsistency between VIIRS and MODIS products [54]. The relationship among VIIRS LAI (x-axis), DLAI (y-axis) and the difference (colorbar) between VIIRS and MODIS LAI retrievals using all composition periods from 2012 to 2016. "DBF" and "ENF" represent "deciduous broadleaf forest" and "evergreen needleleaf forest", respectively.

Limitation and Future Direction
This study investigates the spatiotemporal consistency between S-NPP VIIRS and EOS MODIS LAI/FPAR products. Additionally, the uncertainty of VIIRS products was evaluated using available ground measurements. However, it should be noted that the quantified uncertainty is limited by an availability of the ground measurements (i.e., spatiotemporal and biome-specific coverage). As shown in Figure 8, we can know that the evaluation of VIIRS product uncertainty only achieved stage 1 ("Product accuracy is assessed from a small set of locations and time periods by comparison with in-situ or other suitable reference data." https://viirsland.gsfc.nasa.gov/Val_overview.html) based on CEOS validation stage hierarchy. Thus, adding more ground measurements from spatiotemporally and biome-specific well-distributed field campaigns are required. To achieve a robust evaluation result, we plan to synthesize in-situ reference data from globally distributed research networks (i.e., FLUXNET, Chinese Ecosystem Research Network (CERN), Terrestrial Ecosystem Research Network (TERN), et al.) by deploying a robust upscaling strategy developed by [55] to improve the spatiotemporal coverage of validation. These networks are not originally intended for validation purposes, but they have been continuously collecting in-situ LAI or FPAR Figure 9. (a) The variation of saturation rate with VIIRS LAI retrievals; (b) The relationship among VIIRS LAI (x-axis), DLAI (y-axis) and the difference (colorbar) between VIIRS and MODIS LAI retrievals using all composition periods from 2012 to 2016. "DBF" and "ENF" represent "deciduous broadleaf forest" and "evergreen needleleaf forest", respectively.

Limitation and Future Direction
This study investigates the spatiotemporal consistency between S-NPP VIIRS and EOS MODIS LAI/FPAR products. Additionally, the uncertainty of VIIRS products was evaluated using available ground measurements. However, it should be noted that the quantified uncertainty is limited by an availability of the ground measurements (i.e., spatiotemporal and biome-specific coverage). As shown in Figure 8, we can know that the evaluation of VIIRS product uncertainty only achieved stage 1 ("Product accuracy is assessed from a small set of locations and time periods by comparison with in-situ or other suitable reference data." https://viirsland.gsfc.nasa.gov/Val_overview.html) based on CEOS validation stage hierarchy. Thus, adding more ground measurements from spatiotemporally and biome-specific well-distributed field campaigns are required. To achieve a robust evaluation result, we plan to synthesize in-situ reference data from globally distributed research networks (i.e., FLUXNET, Chinese Ecosystem Research Network (CERN), Terrestrial Ecosystem Research Network (TERN), et al.) by deploying a robust upscaling strategy developed by [55] to improve the spatiotemporal coverage of validation. These networks are not originally intended for validation purposes, but they have been continuously collecting in-situ LAI or FPAR data. These ground data have not been frequently utilized for assessing moderate resolution remote sensing products due to lack of spatial representativeness arising from the scale-mismatch problem. However, a new approach called Grading and Upscaling of Ground Measurements (GUGM) is able to resolve this scale mismatch issue and maximize the use of time-series of ground measurements from these networks. As the approach enables to provide frequent time-series of ground measurements, it would also provide a unique opportunity to evaluate the temporal performance of the LAI/FPAR products [56], allowing better understanding of the structure of the uncertainties and their evolution along different seasonal or annual contexts [33,57]. Therefore, more ground measurements covering the VIIRS product time-period (2012-present), will be collected from the global network of sites and used to better quantify the uncertainties, particularly the performance of time-series attached to the VIIRS product in a future study. This direction will help us to achieve a higher validation stage for VIIRS LAI/FPAR product.
In this study, we focused on evaluating continuity and consistency between retrievals from VIIRS and MODIS, and on quantifying uncertainty via direct validation approach. In general, "indirect validation" comparing VIIRS LAI/FPAR products with products from other moderate resolution sensors can afford one to identify the respective strengths and weaknesses of the products, thus leading to improvements in the next version of products, and to improve the confidence in research-quality of VIIRS products. Moreover, the intercomparison with other similar products is also necessary to achieve stage 2 and stage 3 defined from CEOS validation stage hierarchy ("Spatial and temporal consistency of the product and consistency with similar products has been evaluated over globally representative locations and time periods") [37]. Therefore, we will further analyze VIIRS products using the intercomparison with other global moderate resolution LAI/FPAR products such as GEOV1 [5], GLASS [6] and so on. A series of metrics and qualitative checks proposed by [36] will be adopted to assess the continuity, consistency and precision between VIIRS and other products.
In addition to the provisional results reported in [11,15], our multi-year evaluation clearly corroborates that the implemented heritage MODIS LAI/FPAR algorithm and a physically proven algorithm for achieving inter-sensor consistency allow the generation of multi-sensor LAI/FPAR ESDRs. Our results also satisfy the requirement of stability implying the potential of generating MODIS-like VIIRS LAI/FPAR ESDRs. In light of our S-NPP VIIRS LAI/FPAR experience and knowledge on developing a multi-sensor product, we are confident that we will be able to generate seamless and continuous VIIRS products from JPSS programs with minimal effort. This future plan will build a foundation of continuing sequential JPSS satellite programs, which include already budgeted JPSS-1 to JPSS-4 extending a more than 30+ years of long-term ESDRs.

Conclusions
It is critical that LAI/FPAR products are generated with high accuracy and precision, but more importantly, they must be produced with consistent algorithms across different sensor platforms in order to maintain a continuous and well-characterized data record. This is especially important for the VIIRS LAI/FPAR product on SNPP, which bridges the gap between NASA's EOS satellites and the next generation JPSS platforms. Thus, to evaluate the spatiotemporal consistency between VIIRS and MODIS LAI/FPAR products, multi-year time series of VIIRS and MODIS LAI/FPAR retrievals have been retrieved and their consistency is evaluated in this study. For the spatiotemporal consistency, the global mean differences between VIIRS and MODIS products at a native temporal composite (i.e., 8-day) or seasonally integrated are less than −0.006 for LAI and −0.002 for FPAR, respectively, indicating a reasonably high consistency between two sensors. These estimates meet the stability requirement for long-term LAI/FPAR ESDRs from multi-sensors suggested by GCOS. Generally the performance of the LAI/FPAR consistency for non-forest biomes is better than forest biomes. Additionally, the evaluation also showed a similar interannual variation and a comparable main algorithm retrieval rate from VIIRS with respect to those from MODIS during the study period (2012-2016), although several sources for inconsistency (e.g., stability of input surface reflectance) were discussed. Twenty-one true LAI, 26 effective LAI and 24 FPAR ground measurements with high reliability were used to validate both VIIRS and MODIS LAI/FPAR products. The results indicated that LAI retrievals from both sensors are closer to true LAI than effective LAI because of the clumping correction in the algorithm. We found a comparable LAI uncertainty/relative uncertainty of VIIRS (RMSE = 0.60/42.2%) with respect to that of MODIS (RMSE = 0.55/39.3%) and both VIIRS (RMSE = 0.10/24.4%) and MODIS (RMSE = 0.11/26.0%) showed an overestimation of FPAR. This validation practice also supports the feasibility of producing consistent retrievals from VIIRS and MODIS. However, it also should be noted that more accurate ground measurements and better representation of different vegetation types in different LAI/FPAR ranges are required to refine the uncertainty evaluation of VIIRS LAI/FPAR product. Overall, this study imbues the confidence that S-NPP VIIRS and EOS MODIS can generate a continuous and well-characterized LAI/FPAR ESDRs. This important achievement plays an important role to bridge the gap between NASA's EOS satellites and the next generation of JPSS platforms.
Supplementary Materials: The following are available online at www.mdpi.com/1999-4907/9/2/73/s1, Figure S1: Global land cover map used in this study, Figure S2: The quality control based on "FparLai_QC" and "FparExtra_QC" layers for VIIRS and MODIS products under study, Table S1: Characteristics of the available ground sites and aggregated LAI_effective/LAI_true/FPAR values over 3 km × 3 km area, Figure S3: Spatial distribution of the BELMANIP-2.1 network and ground measurements for VIIRS LAI/FPAR product assessment, Table S2: The mean and standard deviation of LAI/FPAR difference between VIIRS and MODIS products (VIIRS-MODIS) at seasonal scale for each biome type from 2012 to 2016, Figure S4: The examples for the comparison between VIIRS and MODIS products at the native spatiotemporal (500-m/8-day) resolution, Figure  S5: The annual variation of global retrieval rate of different algorithm paths for eight biome types (B1-B8 are same as in Figure S1), Figure S6: The annual variation of mean difference of retrieval index (DRI, calculated from VIIRS-MODIS) for non-forest and forest biomes in January and July from 2012 to 2016.