Revisit the Performance of MODIS and VIIRS Leaf Area Index Products from the Perspective of Time-Series Stability

As an essential vegetation structural parameter, leaf area index (LAI) is involved in many critical biochemical processes, such as photosynthesis, respiration, and precipitation interception. The MODerate resolution Imaging Spectroradiometer (MODIS) and Visible Infrared Imager Radiometer Suite (VIIRS) LAI sequence products have long supported various global climate, biogeochemistry, and energy flux research. These applications all rely on the accuracy of the product's long time series. However, uncontrolled interferences (e.g., adverse observation conditions and sensor uncertainties) potentially introduce substantial uncertainties to time series in product applications. As one of the most sensitive areas in response to global climate change, the Tibet Plateau (TP) has been treated as a crucial testing ground for thousands of studies on vegetation. To ensure the credibility of the studies arising from MODIS/VIIRS LAI products, the temporal quality uncertainties of data need to be clarified. This article proposed a method to revisit the temporal stability of the MODIS (MOD and MYD) and VIIRS (VNP) LAI in the TP, expecting to provide useful information for better accounting for the uncertainties in this area. Results show that the MODIS and VIIRS LAI were relatively stable in time series and available to be used continuously, among which the temporal quality of the MODIS LAI was the most stable. Moreover, the MODIS and VIIRS LAI products performed similarly in both time-series stability and time-series anomaly distribution, magnitudes and fluctuations. The time-series stability evaluation strategy applied to the MODIS and VIIRS LAI can also be employed to other remote sensing products.


I. INTRODUCTION
A S AN essential climate variable recognized by the Global Climate Observing System (GCOS) of the United Nations, leaf area index (LAI) is a basic parameter for measuring the vegetation canopies [1]. LAI is generally defined as the one-sided green leaf area per unit ground horizontal area for broadleaf canopies and as the projected needle leaf area for coniferous canopies [1], [2], [3]. This essential variable plays a key role in hydrology, biogeochemistry, and ecosystem models that connect vegetation to the climate observing system through the carbon, water cycles, and radiation [4], [5], [6], [7]. Retrieving on a clear theoretical basis [8], MODerate resolution Imaging Spectroradiometers (MODIS) LAI products are one of the most long-term LAI products and have been widely used in various fields during the past two decades [9]. MODIS LAI products have attracted the interest of the scientific community due to their successful applications in estimating corn yield [10], simulating terrestrial carbon [11], tracking canopy recovery rates and trajectory after fire [12], analyzing large seasonal swings in leaf area of Amazon rainforests [13], as well as investigating greening or browning of terrestrial vegetation and its drivers during these decades [14], [15]. However, both Terra and Aqua MODIS have far exceeded their design life, they could be terminated in the near future [16], [17]. As the successor of MODIS, the Visible Infrared Imager Radiometer Suite (VIIRS) instrument onboard the Suomi National Polar-orbiting Partnership and Joint Polar Satellite System is expected to continue the long-term Earth system data records (ESDRs) of the MODIS [18].
Considering the dependence of applications on product quality, quality assessment of the product is warranted. Product validation methods can be categorized into three approaches, direct validating using filed measurements, empirical validating with related parameters, such as climatic variables, and intercomparing with other products [19], [20], [21], [22], [23], [24]. Among them, intercomparison can effectively analyze the spatiotemporal consistency of products by comparing differences among existing products [25], generally for a specific period. In this approach, most of the validations usually focus more on the spatial distribution than the temporal variation of the data. As a fact, products with long periods tend to be more widely used, especially for land surface models [26], [27] and the global dynamic monitoring of terrestrial ecosystems [28]. However, This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the long-term datasets were inevitably subjected to undesirable effects due to atmosphere disturbance, sun-sensor geometries, etc. [29]. These effects lead to a lack of temporal consistency of time series, which may introduce considerable uncertainties in the time-series analysis. To obtain reliable results, it is, therefore, necessary to perform the temporal consistency evaluation over the datasets. Existing studies have been conducted to evaluate the temporal consistency of extant commonly used datasets. Tian et al. [30] analyzed the temporal consistency of multisensor Normalized Difference Vegetation Index (NDVI) time series and picked out the GIMMS3g as the most appropriate choice for NDVI trend analysis. Loew et al. [31] investigated the long-term consistency of a surface albedo product and proposed a method to compensate for uncertainties. Brewin et al. (Sea-WiFS) evaluated the temporal consistency of chlorophyll-a in sea-viewing wide field-of-view sensor, MODIS-Aqua and medium resolution imaging spectrometer, concluding that the combining of the three sensors for trend analysis was supported [32].
Several rounds of empirical validations have been conducted to confirm the validity and reliability of the MODIS LAI products. Weiss et al. [33] compared the CYCLOPES LAI with the MODIS LAI, concluding that they had generally consistent seasonality. Fang et al. [20] provided an intercomparison of the MODIS, GEOV1, GLASS, GLOBMAP, and JRC-TIP LAI products between 2003 and 2010, which helps the scientific community better understand their characterizations and uncertainties. Yan et al. [23] assessed the MOD15A2H Collection 6 (C6) product through intercomparisons with GLASS, CYCLOPES, and GEOV1 LAI products, defining the "annual missing data rate" to conduct a temporal comparison over the four products. Meanwhile, the quality of VIIRS LAI products was receiving more and more attention and always being compared with that of the MODIS LAI. Xu et al. [16] analyzed the MODIS and VIIRS LAI products for Spatio-temporal consistency and uncertainty from 2012 to 2016, indicating that the consistency between the two sensors had met the stability requirement for long-term LAI ESDRs from multisensors suggested by GCOS. Yan et al. [34] performed stability analysis on MODIS and VIIRS LAI product quality, concluding that the performance stability of their retrieval algorithm was sufficient to support the trendrelated studies published so far. Nevertheless, the current quality validations of MODIS/VIIRS LAI products have focused more on a given quality flag while ignoring temporal consistency. So far, no study has been made on the evaluation of time-series stability (TSS) for MODIS and VIIRS LAI products. To fill this gap, two TSS indicators, namely TSS and time-series anomaly (TSA), are proposed in this article to revisit the performance of the MODIS and VIIRS LAI from the perspective of TSS.
The TP provides considerable feedback to global environmental changes, which is crucial to reveal environmental change processes and mechanisms on the TP, as well as their influences on and responses to global climate changes [35], [36]. Because of extremely sensitive to climate change, TP is regarded as a crucial testing ground for how humans and the environment collide in a globally warmed world [37]. During these several decades, thousands of studies connected to climate change were conducted in TP. There have been over 1300 vegetation-connected papers about TP published up to 2021 according to the Web of Science. Therefore, dedicating further attention to MODIS and VIIRS LAI products quality from a perspective of time series is extremely significant. For these reasons, we choose TP as the study area, expecting to provide some references for data uncertainties in this area.
Consequently, this article proposed a method for objectively revisiting the TSS of the MODIS and VIIRS LAI as the following structure: Section II describes materials and methods for the TSS revisiting of the MODIS and VIIRS LAI products, as well as the general geography of the study area. Section III details the results about the spatial distribution of LAI and TSS/TSA, the trend and interannual stability of TSS/TSA, fluctuation of TSS within a year and the correlation between TSS and atmospheric conditions/retrieval index (RI). We discuss the correlation between TSS/TSA and LAI retrieval algorithm, the reason for seasonality of TSS, and the overall evaluation of the time-series quality of the MODIS and VIIRS LAI products in Section IV. Finally, Section V concludes this article.  [16], i.e., the VNP15A2H was obtained at almost the same time as the MYD15A2H. The MODIS LAI was retrieved from its main algorithm [based on Radiative Transfer Model (RTM)], or a backup algorithm (based on an empirical relationship between the NDVI and LAI) [2], [38]. As the successor of MODIS, the VIIRS LAI retrieving adopted the same algorithms as the MODIS LAI [38]. Datasets used here all have a spatial resolution of 500 m and a temporal resolution of 8 days. MOD15A2H for 2000, MYD15A2H for 2002, and VIIRS data for 2012 were excluded because they were not complete for those years. Generally, there are 46 composites per year in selected datasets, but owing to the influence of the sensor or other factors, one or two composites may be lost. We have sorted out the lost composites and listed them in Table I attached to Appendix. In this study, we have made up the lost data according to the mean values of other years.

II. MATERIALS AND METHODS
The LAI data were projected on a sinusoidal grid and distributed as standard Hierarchical Data Format files. Each file includes six science datasets (SDSs): Fraction of Photosynthetically Active Radiation absorbed by vegetation (FPAR),  [42] LAI, FparLai_QC, FparExtra_QC, FparStdDev, and LaiStd-Dev. LAI and LaiStdDev layers store the LAI retrieval and its standard deviation. FparLai_QC and FparExtra_QC layers provide information about the algorithm path and atmosphere condition, respectively [39], [40], [41], [42]. In this study, we used the LAI, FparLAI_QC, and FparExtra_QC SDSs. LAI is for calculating TSS, TSA, as well as taken as a reference for analysis. FparLAI_QC is for extracting the retrieval algorithm path and calculating the RI. FparExtra_QC is for extracting cloud and aerosol. We noted that although the SDSs of MODIS and VIIRS are the same in name and storage format, there are still differences in QC files. The values of FparExtra_QC in MODIS and VIIRS LAI are listed below in Tables II and III in the Appendix.
Although there are differences in the aerosols and clouds detecting methods between MODIS and VIIRS, we have tried our best to narrow the differences and keep them at the same level according to the description in user's guides. The specific strategy was detailed in Appendix B.
2) Land Cover Data: The vegetation classification map is the only auxiliary data used in this study. The MODIS land cover product (MCD12Q1) is created by classifying spectral-temporal features derived from reflectance data. It maps the global land cover type at 500 m spatial resolution at an annual time step [43]. Here, we used the LAI legacy classification schemes in which global vegetation is stratified into eight biomes [B1: Grasslands; B2: Shrublands; B3: Broadleaf Croplands (BC); B4: Savannas; B5: evergreen broadleaf forests (EBF); B6: deciduous broadleaf forests (DBF); B7: Evergreen Needleleaf Forests (ENF) [44]; and B8: Deciduous Needleleaf Forests (DNF)] [23], [43]. Considering there is hardly any DNF (B8) occupying the surface of the study area, the following experiment analyzed only the former seven land cover types.
Nineteen years (2001-2019) of MCD12Q1 data were used in this study to control for land cover changes and corresponding effects on LAI magnitude and product quality, thus eliminating the interference brought by land cover changes. We first located the pixels with no land cover change during 2001-2019 and extracted these pixels as mask. According to the mask, we extracted the LAI data with no land cover change during the study period, then examined the extracted data and made up the deficient values with the mean values of other years. The following experiments were conducted with the LAI dataset that has been processed as mentioned.

B. Study Area
The Tibet Plateau (TP) is located in the southwestern of China, covering the extent about 73°19 -104°47 E, 26°00 -39°47 N. Referred to the third pole by scientists, the TP covers approximate five million km 2 with an average elevation of above 4000 m, more than 100 000 km 2 of glaciers included [35], [37]. The land cover type of grasslands spreads the most broadly, which almost occupies half of the surface of the TP [45]. Recent studies have found that grasslands vegetation may contribute considerably to the global carbon cycle [46], [47], [48]. In China, the TP demonstrated a peak ratio of the carbon gap flux to observed net primary productivity, which means strong potential for carbon sequestration [49]. This is certainly significant at a time when the world is striving for carbon neutrality. Besides, there are several forests along the southeastern. The spatial distribution of vegetation cover types in the TP shows the law of zonality. From southeast to northwest, the vegetation cover type changes from forest to pelouse, then to grasslands, and finally to the desert with horizontal zonal distribution [50]. Because of complex terrain and geologic geomorphic structure, the climate in the TP is diverse and complex. The year-round temperature and precipitation of the overall area are relatively low, which differ from parts and distribute not evenly. Contrasting seasons of rain and drought characterize the climate, typically with more precipitation in May to September [51], [52]. Affected by altitude and hydrothermal conditions, the climate type transfers from the warm and wet one in the southeastern to the cold and dry one in the northwestern [53].
We collected vector edge data of the TP from Resource and Environment Science and Data Center (https: //www.resdc.cn/), collected temperature and precipitation data from ERA5 (https://www.ecmwf.int/en/forecasts/datasets/ reanalysis-datasets/era5/) to draw the general situation of vegetation cover types, perennial mean temperature distribution and perennial mean accumulated precipitation in the TP (shown in Fig. 1).

1) TSS Analysis:
In this article, we proposed a metric namely TSS to quantify the fluctuation of time-series. TSS is defined as the distance from the value at the target moment to the linear interpolation line, and the linear interpolation line can be calculated by the data at the former and latter time series of the target moment date.
In the ideal situation, LAI shows the characteristic of monotone increasing and decreasing as the season change. The ideal curve of LAI within a year can be expressed as the smooth green bell curve, as shown in Fig. 2(a). For clearer visual effect and better illustration, we zoom in the part in red rectangular of Fig. 2(a)-(c). Corresponding to the temporal resolution of the experiment data (MOD15A2H/MYD15A2H/VNP15A2H), we extracted discrete LAI values [yellow dots in Fig. 2(a) and (b)] from the ideal curve at an interval of eight days. Therefore, the yellow dots can represent the ideal LAI values with good TSS. However, in reality, the data we can only get is that with jitter due to observation conditions or sensor uncertainties and so on. So assumed that the dark red dots in Fig. 2(c) represent the observation LAI values (i.e., MOD LAI, MYD LAI, and VNP LAI in this article). The dark red dots jumping up and down around the curve denotes the observation LAI with jitter in practice. Linear interpolation lines [bright red dotted line in Fig. 2(b) and (c)] were established between the 16-day-interval data. After that the distance from the data at the middle moment to the line can be calculated, thus gaining absolute TSS of that moment [blue dotted line in Fig. 2(c)]. Ideally, LAI increases or decreases smoothly with no sudden jitter, which leads to the result that the absolute TSS tends to be very small. The phenomenon that the blue dotted lines can be barely seen in Fig. 2(b) illustrates the fact. This is consistent with the law of vegetation growth in the natural state, regarding the ideal LAI data in the curve demonstrating good TSS. Whereas, under the circumstance of data with uncertainties, the jitter will lead to relatively high absolute TSS as the blue dotted lines with certain lengths we can see in Fig. 2(c), which means relatively instability time-series to data. Therefore, the metric TSS we proposed in this article is qualified to depict the jitter of data, i.e., the stability of time-series. The specific calculation process is as follows.
Step 1 Calculating the linear interpolation line: Assumed that the target data X to calculate absolute TSS was gained at the moment t, the two adjacent time-series data were obtained at the moment t −1 (the previous moment) and t 1 (the latter moment), respectively. According to the linear interpolation formula, the linear interpolation line y(t) between the data at the moment t −1 and t 1 should be In which, the X(t −1 ) should be X(i, j, t −1 ) and the X(t 1 ) should be X(i, j, t 1 ) in time-series, because calculation is conducted pixel by pixel, we omit the (i, j) of X(i, j, t) in equations out of brevity and reducing confusion.
Step 2 Calculating the distance from the value at the target moment to the interpolation line: The distance to be calculated, i.e., absolute TSS, can be calculated with the formula for the distance from a point to a line. Substituting the value X(t) at the target moment into the distance calculation formula, absolute TSS can be gained as To figure out how much the data jitter accumulated within a year, we added up all the absolute TSS within a year [i.e., added up all the blue dotted lines in Fig. 2(c)], and then divided by the year number of the data to gain a Multi-year Averaged Yearly Accumulated (MAYA) value, such that In which, m is the composite number within a year and values from 1 to b (for LAI in this article, b is 46); a is the year number, which is 19 in MOD15A2H, 17 in MYD15A2H and 7 in VNP15A2H respectively in this article. The red dotted line separated the data into two parts: the data whose absolute SA ranked the top 90% was put in Part A, while the data whose absolute SA ranked the final 10% of the total data was set in Part B. The data in Part B is defined as outliers. The pixel number of the outliers within a period is TSA.
Step 3 Dividing the absolute TSS by LAI to obtain the Relative TSS: The absolute TSS gained at step 2 represents an original distance from the target moment data to the interpolation line, which is calculated by LAI. Whereas, LAI differs in different land cover types. Besides showing relatively high values in areas with not-so-good time-series data, absolute TSS could also show higher values in the places with higher LAI values.
To eliminate the disturbance of LAI values themselves, we process the absolute TSS by dividing absolute TSS by LAI to get relative TSS. Relative TSS is a unitless indicator that can quantify the relative temporal stability of the LAI among different vegetation cover types, making the time-series quality of different vegetation cover types comparable at the same scale. relative The same as the above, to quantify the data accumulation jitter of the whole year, we calculated the MAYA relative TSS as where the variable meanings are the same as (3).
Both in absolute TSS and relative TSS, high values imply a great deal of data jitter away from the ideal situation, i.e., poor TSS of data; low values imply a stable time-series performance of data conversely.
2) TSA Analysis: In this article, we used the standardized anomaly (SA) to quantify the deviation of LAI in a certain composition from the multiyear average value [38]. SA was calculated as where X(t) and t are the variable time-series and the time (values from 1 to 46 in this article), respectively, mean_t(X(t)) and Std_t(X(t)) meaning the average and standard deviation metrics calculated from the t dimension.
In this article, we counted the distribution of the SA absolute values over the years that the three products cover and drew a histogram (see Fig. 3). According to the histogram, we extracted the pixels whose absolute values of SA exceed the other about 90% pixels (i.e., above 1.65 in this study) as outliers. TSA is so defined as the number of pixels that are beyond the threshold value in the region per year, which quantifies the abnormal values The same as the above, to quantify the anomaly pixels of a year, we summed the TSA in every pixel up within a year, and then divided by the year number to get a MAYA TSA as where the variable meanings are the same as the (3).

3) Trend Analysis:
We used the unary linear regression trend estimation method to recognize the trend of TSS and TSA in this study [54]. The unary linear regression equation can be described as where Data is TSS or TSA, i is the year from 1 to n, and θ slope represents the trend of the data. When θ slope is greater than zero, the data is going on an uptrend ward with the increase of time, and vice versa. The absolute value of θ slope multiplied by ten indicates the rate of the increase or decrease. The greater absolute value of θ slope represents the more acute variation trend. Meanwhile, we employed the Mann Kendall (MK) test for all trend analyses to make sure that the detected trend was statistically significant. The MK test is a common statistical tool for climate diagnostics and prediction [see (10)- (12)], which determines whether there are significant trends in a time-series [34], [55] where S is the sum of the step function values of the difference between the values at any two different points in the time series and n is the number of data, m is the number of tied groups (recurring data sets) in the data, and t i is the number of the ties (the number of repeats in the extent i). When |Z s | > Z 1−α/2 , the null hypothesis (i.e., no trend) is rejected and the α is a special significance level. In this article, we selected α = 0.05 and the Z 1−α/2 = 1.96.

A. Spatial Distribution of LAI and TSS/TSA
To compare the spatial distribution and find the similarity of the LAI, absolute TSS, relative TSS, and TSA, we drew their multiyear averaged values in Fig. 4. Compared with the multicomposite averaged values, multiyear averaged values work in both summarizing the characteristics of each biome and eliminating the influence of the extreme values within a year. Comparing side-by-side [see Fig. 4(1)-(3)], these three products demonstrate similar spatial distribution patterns and share almost the same maximum, minimum and mean values. This is a positive result, especially for VIIRS. As the successor of MODIS, the original aim of VIIRS is to extend the valuable Earth observation dataset of MODIS.  (9), here the percentage in the color bar represents the change rate per decade, numerically equivalent to the slope times 10, then divided by the original value, lastly multiplied one hundred percent. The white in study area represents the data with no trends or no significant trend.
The distributions of LAI and absolute TSS are similar, which both stretch from rather high values in the southeast to rather low values in the northwest. The grasslands perform fairly well in absolute TSS while the forests show high values. But this could have resulted from the LAI values themselves, as high LAI over forests and low LAI over herbaceous vegetation. We can find that the southwestern forests do not show very high values anymore after dividing by LAI to get relative TSS [see Fig. 4(d-1)-(d-3)]. The eastern pelouse shows good TSS. The high relative TSS is mainly concentrated in forests and pelouse, particularly forests. To make the data under different vegetation cover types comparable, we discuss the relative TSS in the following parts only.
As for TSA, most of the parts showed similar distribution. The averaged TSA during a year in each pixel of these products was 3. It means that there were about three anomaly data showing up in each pixel within a year which contains 46 composites. Conspicuously, anomaly pixels in pelouse were detected more in the MOD15A2H [see Fig. 4(d-1)] than in MYD15A2H [see Fig. 4(d-2)], and least in the VNP15A2H [see Fig. 4(d-3)]. Another, the VNP15A2H detected more anomaly pixels in the southwestern forests than MOD15A2H and MYD15A2H.

B. Trend Analysis and Interannual Stability Analysis of TSS/TSA
TSS can be judged from interannual variation and intra-annual variation. Trend analysis and interannual variation of TSS and TSA are presented here. Intra-annual variation of the TSS is discussed in Section III-C.
We extracted the pixels that have passed through the MK test to test the trend. Fig. 5 shows the trend distribution of the relative TSS and TSA. The pixels with increasing trend were marked in green to blue with the increment of increasing degree, while the decreasing ones were marked in yellow to red. Overall, points with trends almost show an uptrend (green or blue, nearly no yellow or red). Some fragmentary points in the northeastern grasslands show an upward trend in relative TSS, which is close to Qilian Mountains. TSA distributes alike. Although there is a sporadic distribution of color spots (points with trends), they are still less compared to the whole quantity and show no region-specific trends, which indicates that the data quality of the three products is relatively stable and with no obvious change trend.
Since the period of VNP15A2H only contains seven years (2013-2019) and this is too short to statistically support trend analysis, we take it as a reference only in this article.
To see more clearly how many pixels increased or decreased in different land cover types over these years, we calculated the proportion of pixels with increasing or decreasing trends in each land cover type. Overall, there are more pixels with a trend in the TSS [see Fig. 6(a)] than the TSA [see Fig. 6(b)], more pixels with an upward trend than a downward trend. The MOD TSS/TSA take up a greater percentage of increasing trend pixels than the MYD TSS/TSA. Approximately 22.6% of MOD TSS pixels of the BC land type show an increasing trend, which makes BC stand out in the bar chart of Fig. 6(a). BC also shows a tendency to increase more in MYD TSS (11.9%), MOD TSA (14.5%), and MYD TSA (12.2%) compare to other land cover types. Although the data for this land cover type shows the least stability compared to other land cover types, it is still relatively stable overall with the rest about 80%TSA pixels and 90% TSA pixels showing no trends. Pixels with an increasing trend of shrublands TSS in MOD and MYD are almost at the same proportion, which is about 10% and slightly higher than others. The percentages of trending pixels in other land cover types are all below 10%. The statistics show that the data remain stable  in the vast majority of land cover types in terms of time series stability.
Likewise, we do not analyze the trend of VNP because of its short data period.
The yearly averaged relative TSS in different vegetation cover types show a characteristic that fluctuates but no trends in most of them. The three products show a rather consistent wave in the same vegetation cover type both in the TSS and the TSA. Because of the high coverage of grasslands, the overall picture displays similar fluctuation with the grasslands. Only relative TSS in Shrublands [see Fig. 7(b)] has passed the MK test and shows a slight upward trend. TSS in the other land cover types have no trends and demonstrate nearly stability, which confirms the continuous availability of the products. On the whole, the broken lines of MOD TSS always lie below the other two, which means that the MOD15A2H temporal stability is better than that of MYD15A2H and VNP15A2H. That is to say, the MOD15A2H product is with little jitter than the other two and closer to the ideal value. The MYD15A2H and VNP15A2H stay at a similar level to TSS. In terms of land cover types, the As we can see from Fig. 7, although the MODIS has served for over 20 years, which is far beyond the design life, it still shows great vitality and stability. Its long time series and stability will continue to feed the scientific community.
Likewise, the fluctuations of the three products match well with each other. It is clearly shown that anomaly data in VNP15A2H is less than both in MOD15A2H and MYD15A2H. Both of the MODIS products have a comparable number of anomalies. To be specific, the mean number of anomaly data over the whole study area are 3.75, 3.69, and 3.20 for MOD15A2H, MYD15A2H, and VNP15A2H, respectively.
Among the seven land cover types, trends are detected in grasslands, shrublands, BC, and DBF. The overall TSA [see Fig. 8(h)] goes with an upward trend, whose shape and growth rate is close to that of grasslands [see Fig. 8(a)] owing to the high coverage of grasslands. The shrublands contain the least TSA while the BC contains the most TSA. Meanwhile, the TSA growth rate of the BC is beyond other land cover types. The slope of the BC TSA fitting line implies that it could be about more 1.8 and 1.7 anomaly data showing up in every pixel for MOD15A2H and MYD15A2H during the next decade. Consequently, although the growing trend were detected in the TSA, the increments were not significant.
TSA tells the number of anomalies within a year. These anomalies could be induced from the deviation of LAI products themselves or climate. To be clear, it was caused by surface conditions. Considering the TP is getting warmer and wetter during the five decades [37], [56], [57], [58], extreme weather due to climate changes may explain the increment of the anomaly. The TSA indicator thus reflects the climate change of the TP from a side angle.

C. TSS Fluctuation Within a Year
To see specifically how the TSS changes within a year, we plot the TSS values day by day using the adjacent 8-day LAI.
As expected, the TSS of MOD15A2H is slightly better than that of MYD15A2H. TSS of VNP15A2H and TSS of MYD15A2H demonstrate close values, but that of the VNP15A2H with some more sharp mutations. It is interesting to find that the TSS fluctuates significantly during the growing season, which indicates the seasonality of data quality. Data quality being more likely to be unstable in the growing season may result from the increased probability of clouds, aerosols, and precipitation in summer, as will be further discussed later. Meanwhile, forest biomes show relatively high TSS values in the summer season over other biomes. The surface reflectance saturation in dense canopies in summer could have also contributed to the time-series instabilities of data [38]. Despite changing dramatically during the growing season, the TSS of shrublands remains the least among others. TSS of grasslands fluctuates relatively smoothly. This rule tells the different sensitivities of LAI to surface reflectances. The saturated part and the unsaturated part will introduce different uncertainties to the LAI data, thus in turn affecting the stability of data quality [25].

D. Correlation Between TSS and Meteoro/RI
To confirm the relationship between the atmospheric variables and the TSS, we conducted correlation analysis for TSS and aerosol/cloud. Added that according to prior knowledge, poor observation conditions (cloudy, humid atmosphere, etc.) could lead to problematic retrieval, we also conducted correlation analysis for TSS and the RI. The RI is the percentage of pixels retrieved by the radiative transfer equation (RTE)-based main algorithm of the pixels retrieved by both the main and backup algorithm, which characterizes the proportion of the high precision retrievals with good quality at a regional scale [16], [23], [34], [38] where N 1 is the number of pixels retrieved by the main algorithm in study area and N 2 is the number of pixels retrieved by the backup algorithm. A ratio approach was conducted to eliminate magnitude effects. We used the percentage of aerosols and clouds to represent the atmospheric conditions. Correlation analysis of aerosol/cloud ratios and relative TSS in the region was then conducted. Product TSS has little to do with aerosols with the correlation between TSS and aerosol neither high nor significant [see Fig. 10(a)]. Significant correlation is detected between relative TSS of all the three products and clouds (R>0.7) [see Fig. 10(b)]. Surely there were no obvious trends of aerosols and clouds during the two decades according to the color points distribution. This result advance illustrates that product time series stability can be influenced by atmospheric conditions. However, considering there were no trends in the atmosphere variables, we believe that atmospheric conditions have a uniform effect on data quality, which explains the lack of trends in TSS for most land cover types. According to the statistical result in Fig. 10(b), MYD15A2H detected more clouds than MOD15A2H, which is also consistent with meteorological statistics that morning time with fewer clouds than afternoon reported by the National Meteorological Information Center (http://data.cma.cn). This result implies that fewer interferences by the clouds and aerosols can explain parts of the reason that the MOD15A2H product showed better TSS.
As expected, the relative TSS is negatively correlated to the RI [see Fig. 10(c)], which reveals that a higher proportion of the main algorithm led to lower relative TSS, i.e., better TSS. Noted that apart from being affected by the atmospheric conditions, the RI is also infected by the land cover types, as to be discussed in the following part. Fig. 10. Correlation between relative TSS and (a) aerosol/ (b) cloud/ (c) RI. The x-axis represents the mean relative TSS at the pixel level; the y-axis in (a) and (b) represents the annual percentage of aerosol/clouds in the region (The aerosols/clouds from each image within the same year were adding up, and divided the result by the total number of pixels during a year to obtain the annual aerosol/cloud percentage in the region). The RI is calculated as (13). The color bar on the right tells about the time of each point. Blue to red represents the early years to recent years.

A. Correlation Between TSS/TSA and LAI Algorithm Path
The essential difference between the TSS and TSA requires further clarification here. TSS quantifies the TSS of products, which depicts volatility characteristics of product quality and is subject to the LAI product quality. For instance, meteorological conditions (e.g., cloud and aerosol) at the time of observations can affect the product quality and, in turn, affect TSS value. The TSA calculates the number of anomalies within a year, which will be affected by the LAI value itself. The anomalies could be induced by product uncertainties, climate or LAI anomaly of the surface itself, etc. So the detected trend of TSA is possibly due to increased LAI anomalies caused by extreme weather due to climate change (see Fig. 8). Based on the above, we discuss several scenarios that may lead to fluctuations in TSS and TSA from the LAI algorithm mechanisms standpoint.
There are three algorithm paths in retrieving MODIS and VIIRS LAI: the main algorithm, a look-up table with no saturation and saturation [39], [59], and the backup algorithm [3], [60]. Their LAI retrieval accuracy decrease in order [22]. The product accuracy will be high with a higher proportion of main algorithms. RI is, therefore, defined as the percentage of pixels for which the main algorithm produces a retrieval, and it characterizes the overall quality at a regional scale. During the growing season, cloudy or high aerosol atmospheric conditions all add to the difficulties of estimating LAI using the main algorithm, which lowers the RI in the study area. Added that the RI is sensitive to high LAI values. To be specific, the domain of the main algorithm retrievals is more limited in dense vegetation with higher near-infrared band reflectance and lower red band reflectance [22]. Due to the interannual variability of surface reflectance, the limited retrievals will introduce fewer good quality and high precision retrievals [16].
On the other hand, the saturation effect will work to lower the retrieval accuracy when the LAI gets higher because the surface reflectance cannot contain enough information [61]. According to the spatial distribution, we found it noticeable that data quality differs from the area. Generally, the performance of the LAI TSS for nonforest biomes is better than forest biomes. We believe that was an outcome of the combination of low RI and saturation effects. They are both unfriendly to high LAI value regions. Either low RI or saturation indicates low precision of products, which could further contribute to the volatility of products. This may explain the reason why most of the data with poor TSS concentrates in the forest region. Added that the misclassification of biome types will also introduce LAI uncertainty.

B. Explanation of TSS's Seasonality
Correlation between the TSS and atmospheric elements (see Fig. 10) talks about the data temporal stability being influenced by weather conditions. The effect of atmospheric conditions on the data TSS happened most during the growing season (see Fig. 9). Because in the growing season (generally referred to summer), the probability of clouds, aerosols and precipitation is higher. Apart from the aerosols and clouds, we have also conducted a correlation analysis between TSS and precipitation/temperature. We found that there is a significant correlation between TSS and precipitation (MOD: R = 0.5969; MYD: R = 0.5047) whereas the temperature is not, which may be resulted from LAI subjection to the precipitation. The correlation between precipitation and LAI had also been reported by other scholars [23], [62]. We do not specify it for concise. Poor atmospheric conditions during the growing season can reduce the RI, thus lowering the TSS of data. Also, along with the increase of LAI during the growing season, the saturation effect will affect the accuracy of the product to some extent. Both of them could lead to the seasonality of TSS.
Despite the correlation between aerosols and clouds, there were no trends in them, corresponding to the homogeneous effect of meteorological conditions on TSS, this will not lead to a trend in TSS and agrees with Yan et al. [34]. The data quality metrics showing no trends demonstrates the stability and continued availability of the MODIS and VIIRS LAI again. Meanwhile, noted that the aerosol and cloud of VNP are distinguished. This is because the MODIS and VIIRS do not share the quality control methods and files as mentioned in Section II. Because of the absence of CO 2 and H 2 O absorption channels, the cloud detection and flagging mask in VIIRS are limited, in turn such discrepancy results in unlike LAI retrievals [63]. In comparison, the MODIS owes far more bands than the VIIRS, which helps a lot in aerosol and cloud detection, so the aerosol and cloud detection in the MODIS are more reliable. Both the Suomi-NPP and Aqua-MODIS are afternoon-overpassing satellite, their transit time is identical [16], so VNP15A2H and MYD15A2H should have shared the alike level of the cloud and aerosol. However, according to the statistical result in Fig. 10(b), MYD15A2H detected more clouds than VNP15A2H. Meanwhile, it is no denying that different observing conditions in the two sensors would also interfere with the difference we observed. For instance, different selection of daily best input observation may cause a dissimilar input surface reflectance [16], which further causes the difference of LAI retrieval. Theoretically, the clouds of MYD15A2H and VNP15A2H should be more than that of MOD15A2H, because Terra-MODIS passes territory in the morning with little cloud proportion. Statistically, morning time with less cloud than afternoon is in accordance with nature (http://data.cma.cn). As the statistical result in Fig. 10(b), the clouds of MYD15A2H were indeed more than that of the MOD15A2H, while the VNP15A2H detected fewer clouds than MOD15A2H. This is because the cloud may be leak-detected in the VNP products due to the absence of CO 2 and H 2 O absorption channels in the VIIRS [63]. Improving the quality and TSS of products under cloudy or high aerosol atmospheric conditions is still challenging [64], [65].

C. Overall Evaluation of the MODIS and VIIRS LAI Quality
According to the spatial distribution, trend and fluctuation of the TSS and TSA of the three products, we found that all of the three products were stable and performed similarly in temporal quality from the aspects mentioned above. Among the three products, MOD15A2H demonstrated the lowest TSS while the VNP15A2H counted the least TSA. Generally, the TSS of the non-forest biomes performed better than that of the forest biomes. The land cover type of BC detected more TSA than other land cover types, with no significant differences in the amount of the TSA among other six biomes. Massive TSA in the BC could be the result of artificial intervention in surface conditions, such as cultivating and reaping. Considering there were more and full validations on MOD15A2H [16], [23], [25], [34], [66]. Added that the MOD15A2H product owes the longest time series over the three products. It is most recommended to use the MOD15A2H product over the overlapping period.

D. Applicability and Limitations
Overall, uncertainties can be divided into two categories, systematic uncertainties and random uncertainties. Systematic uncertainties are usually caused by sensors or calibration, while random uncertainties depict uncertainties usually introduced by retrieval algorithm. Stability defined by GCOS refers to 'the maximum acceptable change in systematic error measured in decadal time scales' [67], which is hard to quantify from the products since the land surface usually cannot be considered as no-change in a decade. With no specific quality requirement for random uncertainties currently, Fang et al. [68] have referred to the GCOS stability requirement as a proxy quality requirement for temporal uncertainties. In TSS and TSA calculations, both of them eliminate the trend component indicated by average LAI, which means TSS and TSA measure random uncertainties rather than systematic uncertainties, and thus, essentially differ from the stability defined by GCOS. However, we can check the decennary systematicness of the short-time random uncertainties by calculating the trend of TSS or TSA.
With regard to temporal scale, GCOS focus on long timeseries interannual stability (decadal time scale), while TSS and TSA are calculated from every 8-day composition. TSS and TSA mainly measure intra-annual short-time data jitters, but can as well as compare the interannual jitter by summing up all the components within a year. Through intra-annual stability, anomalous jitters introduced by algorithm or atmospheric conditions can be checked and quantified. Note that there should be no long-term trend (systematic uncertainties) of annual TSS or TSA if the sensor and the calibration performance are stable. In another words, the long-term trend of TSS and TSA can be used to quantify the systematicness of product quality.
Admittedly there are some unavoidable limitations in this article. On one hand, as it is known to us all that air in the TP is rarefied, which weakens aerosol scattering and makes the statistics of aerosol fail to be highly representative for other areas. On the other hand, apart from lacking DNF samples, the number of pixels occupying by BC (no more than three thousand) is far less than other biomes. TSS of data in both of the two biomes requires further investigation.

V. CONCLUSION
Data quality underpins the studies. Well-consistency data selection is the first step toward successful research. Data with poor time series stability may even lead to erroneous trends and results, especially in research conducted over long time series. The reliability of studies based on the MODIS and VIIRS LAI depends on their TSS, for which detailed research has been missing to date.
To accomplish this assessment, this study proposed a method to objectively revisit the TSS of MODIS and VIIRS LAI in the TP. We compared the spatial distribution patterns of LAI and TSS/TSA, analyzed the trend and interannual stability of TSS/TSA, observed the TSS fluctuation within a year and investigated the correlation between TSS and the atmospheric elements/RI. The results showed that the TSS of both the MODIS and VIIRS LAI had no trends overall and were relatively stable, which confirmed the continuous usability of the products. Product TSS varied seasonally because of the low RI and saturation caused by poor atmospheric conditions during the growing season. Moreover, results also indicate that although clouds can explain a large part of the poor TSS, the interannual variations in atmospheric conditions with no trends will not result in large changes in product qualities. As for TSA, increasing trends were detected in some biome types, which implied that climate change has slow but notable impact on the surface in the TP. Last, although we not aimed at evaluating the consistency between the MODIS and VIIRS LAI, we found that the MODIS and VIIRS LAI products performed similarly in both TSS and TSA distribution, magnitudes and fluctuations. Temporal quality varying similarly both in spatial and temporal, we can know that the TSS of the MODIS and VIIRS LAI was comparable.
In addition to filling the research gap in the time series evaluation of MODIS and VIIRS LAI data, this study provides a feasible route for selecting data with better time stability. It also can be applied to assess the TSS of other remote sensing products. Besides, our study offered a seminal idea that quality trends can also be used to capture anomalies of the ground surface.

A. Lost Composites of MODIS LAI
According to statistics, there are four composites lost during 2001 and 2019 in MOD15A2H products while there are no lost composites in MYD15A2H and VNP15A2H. The lost composites of MOD15A2H in the study area are listed below in Table I.

B. Strategy of Extracting Aerosols and Clouds From MODIS and VIIRS LAI
According to the user's guide of MODIS and VIIRS LAI, we extracted the aerosol and cloud in MODIS LAI when the No.3 and No.5 bit-string is 1, while we extracted the aerosol and cloud in VIIRS LAI when the No.0-1 bit-string is 11 and No.4-5 bit-string are 10 and 11 (both marked in bold in Tables  II and III). Noted that we assumed that the aerosol in VIIRS LAI marked average or high levels (No.4-5 bit string are 10 and 11) correspond to the detected aerosol in MODIS LAI (No.3 bit string is 1), because the illustration in MODIS LAI tells the detected aerosols should be in average or high level (marked in red in Table II).