Inﬂuence of Snow on the Magnitude and Seasonal Variation of the Clumping Index Retrieved from MODIS BRDF Products

: The foliage Clumping Index (CI) is an important vegetation structure parameter that allows for the accurate separation of sunlit and shaded leaves in a canopy. The CI and its seasonality are critical for global Leaf Area Index (LAI) estimating and ecological modelling. However, the cover of snow tends to reduce the reﬂectance anisotropy of the vegetation canopy and thus probably inﬂuences CI estimates. In this paper, we investigate the inﬂuence of snow on the magnitude and seasonal variation of the CI retrieved from Moderate-resolution Imaging Spectroradiometer (MODIS) Bidirectional Reﬂectance Distribution Function (BRDF) products based on ﬁeld-measured CI and statistics from the global MODIS CI product. We ﬁnd that the backup algorithm can effectively correct abnormally large CI values and obtain more reasonable CI retrievals than the main algorithm without any constraints in snow-covered areas. Validation indicates that the time-series CI product shows the potential in investigating the trajectories of the clumping effect in snow seasons. For evergreen forests, the clumping effect is relatively stable throughout the year; however, for deciduous vegetation types, CI values tend to display signiﬁcant seasonal variations. This study suggests that the latest version of the global MODIS CI product, in which the backup algorithm is used to process the snow-covered pixels, has improved accuracy for CI retrievals in snow-covered areas and thus is probably more suitable as the input parameter for ecological and meteorological models.

The relationship between the Normalized Difference between Hotspot and Darkspot (NDHD) and the CI, which was originally developed by Chen et al. [4], has been used to produce regional and global CI products based on various multi-angular data from the Polarization and Directionality of the Earth's Reflectances (POLDER,~6 km resolution) sensor [4,[32][33][34], Moderate Resolution Imaging Spectroradiometer (MODIS,~500 m resolution) sensor [24,26,[35][36][37][38][39] and Multi-angle Imaging SpectroRadiometer (MISR,~275 m resolution) sensor [5]. The early global CI maps were usually generated with only one annual group of CI values that were retrieved by extracting either the minimum [4,24] or median CI [26]. The CI is usually assumed to be seasonally invariant in land surface models [40][41][42] because of the limited data availability. However, the seasonality of the CI is critical for LAI estimating and ecological modelling [37]. Several studies based on the ground measurements of long-term time series CI data indicate that the CI could have very strong seasonal variations due to species phenology [5,23,43,44]. In addition, various studies based on spaceborne CI products indicate that the retrieved CI values exhibit distinguishable seasonal variations [4][5][6][37][38][39]. Therefore, the seasonality of the CI should be carefully explored in spaceborne CI products.
One difficulty in capturing the seasonality of spaceborne CI products is that the seasonal variation of the CI may be contaminated by short-term fluctuations in the retrieved CI time series data [37]. For example, abnormally small or large values can be found in the CI retrievals from MODIS BRDF products at validation sites [37,38]. Recent studies found that outliers might be caused by ephemeral rain [37] or snow events [39]. Jiao et al. [39] reported that ephemeral snow during winter decreases the accuracy of MODIS Bidirectional Reflectance Distribution Function (BRDF) retrievals during leaf-off seasons, which leads to lower CI quality. Therefore, the influence of snow on spaceborne CI products should be further investigated to improve the accuracy of the global MODIS CI product when investigating the seasonality of the CI.
In this study, we mainly aim to investigate (1) the influence of ephemeral snow on the magnitude of the CI and (2) the seasonal variation of the CI after correcting for or excluding the influence of ephemeral snow. For these purposes, we collect ground measurements of CI time series data in three validation sites and the MODIS snow cover product to analyze the influence of snow on the CI product retrieved by Jiao et al. [39] from MODIS BRDF products. In addition, we summarize and analyze the CI statistics for various International Geosphere-Biosphere Programme (IGBP) classes after correcting or excluding the snow-covered pixels on the global scale to investigate the influence of snow on the magnitude and seasonal variation of spaceborne CI products.

Global MODIS CI Product
The global MODIS CI product used in this study, hereby termed CI MODIS , is retrieved by Jiao et al. [39] with a spatial resolution of~500 m and an 8-day temporal resolution from January 2001 to December 2013. CI MODIS is retrieved from the MODIS BRDF products [45] based on the following linear relationship between CI and NDHD: where A and B are the linear regression coefficients determined by a set of model simulations based on the 4-scale geometrical optical model [46] in Chen et al. [4]. The values of A and B vary with solar zenith angle (θ), spectral band and crown shape [4]. NDHD is defined as follows: where R Hot and R Dark are the reflectance in hotspots and darkspots, respectively. NDHD can be calculated based on the three BRDF parameters in the MODIS BRDF product using the hotspot-corrected version of the kernel-driven model (RTCLSR) developed by Jiao et al. [47]. Here, the MODIS BRDF parameters in the red band are selected to calculate the NDHD since He et al. [26] reported that the CIs retrieved from the red band were more accurate than those retrieved from the near-infrared (NIR) band. Accordingly, CI can be calculated based on Equation (1) using the corresponding coefficients (when θ = 45 • ) in the red band. The abovementioned algorithm based on the linear CI-NDHD relationship is called the main algorithm. If CI retrievals by the main algorithm are outside of the closed interval [0.33, 1.00] or the pixels are covered by snow, then a backup algorithm will be used instead of the main algorithm to achieve more reasonable CI values [39].
The backup algorithm is developed based on the angular index Anisotropic Flat Index (AFX) [48][49][50][51]. The AFX indicates the general variability of BRDF shapes (dome or bowl-shaped BRDFs) in terms of the value around unity. Jiao et al. [39] assumed that the pixels with similar canopy structures would capture similar BRDF shapes that can be identified by similar AFX values (i.e., in an AFX interval of 0.01). Therefore, a lookup table (LUT) between the CI and AFX for each IGBP class can be built based on the accumulated high quality MODIS BRDF data, which are marked as best quality by the MCD43A2 product and are not covered by snow. In the LUT, the CIs for each corresponding AFX value in each IGBP class can be calculated based on the following procedure. Firstly, a density separation method is adapted to slice the AFX range (i.e., 0.5 to 1.5) at an interval of 0.01. Then, the high-quality BRDF shapes in each AFX interval are averaged to achieve a typical BRDF archetype (i.e., an average BRDF shape) that can represent the BRDF shapes in the corresponding AFX interval. Finally, the CIs for each corresponding AFX value in each IGBP class can be calculated from the BRDF archetype using the liner relationship between the CI and NDHD. More details can be found in Jiao et al. [39].
CI retrievals by the main algorithm are marked with a MODIS CI QA flag of 0 or 1, and those by the backup algorithm are marked with a MODIS CI QA flag of 2.

MODIS Snow Cover Product
The MODIS Snow Cover product (MOD10A2, version 6) [52] is derived from the Normalized Difference Snow Index (NDSI). It reports the maximum snow cover extent during an eight-day period in 1200 km × 1200 km tiles with a spatial resolution of~500 m. The MOD10A2 product from January 2001 to December 2013 and other MODIS products used in this study were downloaded from https://search.earthdata.nasa.gov/.
In this study, the dataset of maximum snow cover extent in the MOD10A2 product is used to select or mask the pixels that are covered by snow. The dataset is generated from the MOD10A1 tiles. If snow is observed in a pixel on any day in the period of MOD10A2, the pixel will be marked as snow. If no snow is found, the pixel will be marked as the land cover that is observed most often (snow-free land, lake and so on). If the pixel is covered by clouds for all eight days in the period, the pixel will be marked as a cloud.

MODIS Land Cover Type Product
The MODIS land cover type product (MCD12Q1) describes the land cover properties derived from the observations spanning a year of MODIS data [53]. It incorporates five land cover classification schemes derived using a supervised decision-tree classification method. The IGBP global vegetation classification dataset is used in this study for consistency with the MODIS CI product.

Collected Ground Measurements of CI Time Series
The ground measurements of the CI time series for three validation sites were collected from previous studies [6,37] and compared with CI MODIS to analyze the seasonal variation of the CI. Details on the three validation sites are shown in Table 1. The TP39 and TP74 sites representing Evergreen Needleleaf Forest (ENF) are in Ontario, Canada. The areas are dominated by mature white pine [54]. The CIs at the two sites are measured by the Tracing Radiation and Architecture of Canopies (TRAC) instrument using the gap size distribution [1,25,55,56] method [37]. The Honghe Farm site is in the Heilongjiang Province, China. It is dominated by large, flat, homogeneous rice paddy fields (homogeneity > 5 km). The paddy rice is grown and transplanted in late May; flowering occurs in early July, grain-filling in early August and maturity in early September. The CI of this site was continuously measured from June (Day of Year (DOY) 163) to September (DOY 261) 2012, which includes the main period of paddy rice development from germination to maturity [57].

Analytical Method Regarding the Influence of Snow on the Maginutde and Seasonal Variation of CI MODIS
First, we compare the ground measurements of the CI time series for the three validation sites using CI MODIS to investigate the influence of snow on the magnitude of CI MODIS . The MOD10A2 product is used to mask snow-covered pixels. The CI MODIS of the sites will be recalculated using the corrected IGBP type if the IGBP class in the MCD12Q1 product is not consistent with the real IGBP class of the sites. Then, we investigate the extent of the areas influenced by snow for various IGBP classes based on the global snow cover and CI maps on a specific day in February 2006. Finally, we compare the temporal sequences of the average snow-free CI MODIS with the snow-covered CI MODIS for the typical IGBP classes in the northern hemisphere at latitudes ≥ 30 • N from January 2001 to December 2013 to further investigate the influence of snow on the magnitude and seasonal variation of CI MODIS . Figure 1 shows the seasonal variation of the in situ CI (red circles) and CI MODIS . Figure 1a,b shows that although CI MODIS retrieved by the main algorithm matches quite well with the field-measured CI at the TP39 and TP74 sites on some of the snow-free days, the CI MODIS fluctuates too much with comparison to the in-situ measurements in snow-free seasons. The CI MODIS values retrieved by the main algorithm are much larger than the field-measured CI values on snow-covered days. In contrast, the CI MODIS values retrieved by the backup algorithm can correct the influence of snow and show good consistency with the field-measured CI on snow-covered days. Several CI MODIS values retrieved by the main algorithm are larger than 0.6 or smaller than 0.5 during the period from May to October. The fluctuation of CI MODIS may be caused by noise from rain events [37]. Interestingly, the backup algorithm shows a higher accuracy than the main algorithm when those outliers (CI > 0.6 or CI < 0.5) are retrieved by the main algorithm. In contrast, the main algorithm performs better than the backup algorithm on the snow-free days expect for the days when outliers are retrieved by the main algorithm.

Influence of Snow on CI MODIS at Validation sites
At the Honghe Farm site, a rapid change in the CI can be observed from the in situ measurements of the CI (Figure 1c), since the structure of paddy rice varies rapidly during the growing season. CI MODIS retrieved by the main algorithm significantly captures the minimum, maximum and trajectories of the field-measured CI. The values of CI MODIS are much smaller than the values of the field-measured CI on DOYs 163, 164, 259, and 260 since CI MODIS is retrieved from the MODIS BRDF product, which assumes that the BRDF of the observed target remains consistent during the 16 days. Therefore, it is difficult for CI MODIS to capture the sudden real-time change of the CI. Apparently, the main algorithm performs better than the backup algorithm on snow-free days, which is consistent with our previous result that the main algorithm outperforms backup algorithm except for CI outliers [39]. In contrast, the CI MODIS retrieved by the main algorithm on snow-covered days are all larger than 1.00 and thus should be considerably larger than the real CI since the values of the CI are unlikely to exceed 1.00 in theory [4]. The results indicate that the CI MODIS retrieved by the backup algorithm should be more reasonable than that retrieved by the main algorithm on snow-covered days. The accuracy of the CI MODIS retrieved by the backup algorithm on snow-covered days needs further evaluation with sufficient validation data on snow-covered days. At the Honghe Farm site, a rapid change in the CI can be observed from the in situ measurements of the CI (Figure 1c), since the structure of paddy rice varies rapidly during the growing season. CIMODIS retrieved by the main algorithm significantly captures the minimum, maximum and trajectories of the field-measured CI. The values of CIMODIS are much smaller than the values of the field-measured CI on DOYs 163, 164, 259, and 260 since CIMODIS is retrieved from the MODIS BRDF product, which assumes that the BRDF of the observed target remains consistent during the 16 days. Therefore, it is difficult for CIMODIS to capture the sudden real-time change of the CI. Apparently, the main algorithm performs better than the backup algorithm on snow-free days, which is consistent with our previous result that the main algorithm outperforms backup algorithm except for CI outliers [39]. In contrast, the CIMODIS retrieved by the main algorithm on snow-covered days are all larger than 1.00 and thus should be considerably larger than the real CI since the values of the CI are unlikely to exceed 1.00 in theory [4]. The results indicate that the CIMODIS retrieved by the backup algorithm should be more reasonable than that retrieved by the main algorithm on snow-covered days. The accuracy of the CIMODIS retrieved by the backup algorithm on snow-covered days needs further evaluation with sufficient validation data on snow-covered days.

Areas Influenced by Snow at a Global Scale
In winter, snowfall can often be observed in temperate and frigid zones. Figure 2 shows that snow covers most areas of the North Frigid Zone and the northern region of the North Temperate Zone on DOY 041, 2006. Statistics indicate that nearly 31 percent of the areas in the northern hemisphere (Latitude > 0°N) on DOY 041 are covered by snow. Detailed statistics for various IGBP classes are shown in Figure 3. The ENF and Deciduous Needleleaf Forest (DNF) are more prone to be affected by snowfall. By contrast, Evergreen Broadleaf Forest (EBF) is less prone to be affected by snowfall, since most EBF are distributed in tropical and subtropical regions (where snowfall rarely occurs). Therefore, the influence of snow on CIMODIS should be carefully investigated in winter, especially for needleleaf forests.

Areas Influenced by Snow at a Global Scale
In winter, snowfall can often be observed in temperate and frigid zones. Figure 2 shows that snow covers most areas of the North Frigid Zone and the northern region of the North Temperate Zone on DOY 041, 2006. Statistics indicate that nearly 31 percent of the areas in the northern hemisphere (Latitude > 0 • N) on DOY 041 are covered by snow. Detailed statistics for various IGBP classes are shown in Figure 3. The ENF and Deciduous Needleleaf Forest (DNF) are more prone to be affected by snowfall. By contrast, Evergreen Broadleaf Forest (EBF) is less prone to be affected by snowfall, since most EBF are distributed in tropical and subtropical regions (where snowfall rarely occurs). Therefore, the influence of snow on CI MODIS should be carefully investigated in winter, especially for needleleaf forests. At the Honghe Farm site, a rapid change in the CI can be observed from the in situ measurements of the CI (Figure 1c), since the structure of paddy rice varies rapidly during the growing season. CIMODIS retrieved by the main algorithm significantly captures the minimum, maximum and trajectories of the field-measured CI. The values of CIMODIS are much smaller than the values of the field-measured CI on DOYs 163, 164, 259, and 260 since CIMODIS is retrieved from the MODIS BRDF product, which assumes that the BRDF of the observed target remains consistent during the 16 days. Therefore, it is difficult for CIMODIS to capture the sudden real-time change of the CI. Apparently, the main algorithm performs better than the backup algorithm on snow-free days, which is consistent with our previous result that the main algorithm outperforms backup algorithm except for CI outliers [39]. In contrast, the CIMODIS retrieved by the main algorithm on snow-covered days are all larger than 1.00 and thus should be considerably larger than the real CI since the values of the CI are unlikely to exceed 1.00 in theory [4]. The results indicate that the CIMODIS retrieved by the backup algorithm should be more reasonable than that retrieved by the main algorithm on snow-covered days. The accuracy of the CIMODIS retrieved by the backup algorithm on snow-covered days needs further evaluation with sufficient validation data on snow-covered days.

Areas Influenced by Snow at a Global Scale
In winter, snowfall can often be observed in temperate and frigid zones. Figure 2 shows that snow covers most areas of the North Frigid Zone and the northern region of the North Temperate Zone on DOY 041, 2006. Statistics indicate that nearly 31 percent of the areas in the northern hemisphere (Latitude > 0°N) on DOY 041 are covered by snow. Detailed statistics for various IGBP classes are shown in Figure 3. The ENF and Deciduous Needleleaf Forest (DNF) are more prone to be affected by snowfall. By contrast, Evergreen Broadleaf Forest (EBF) is less prone to be affected by snowfall, since most EBF are distributed in tropical and subtropical regions (where snowfall rarely occurs). Therefore, the influence of snow on CIMODIS should be carefully investigated in winter, especially for needleleaf forests.    Table 2.
As shown in Figure 2, most CIMODIS values for snow-covered areas are either larger than 1.00 or have no value. Statistics indicate that in areas with latitudes ≥ 30°N in the northern hemisphere, 94% of the pixels with a CIMODIS ≥ 1.00 are covered by snow. The snow on the ground is highly isotropic scattered, thus reducing the contrast between the sunlit and shadowed components of the vegetation canopy [58][59][60], resulting in a smaller NDHD. Therefore, a larger CIMODIS will be retrieved from the snow-covered pixels. The results further demonstrate that snow cover reduces the quality of the retrieved CIMODIS, resulting in either abnormally large CI values or no values. Therefore, the CIMODIS retrieved from snow-covered pixels should be corrected or excluded for various applications when studying the seasonal variation of CIMODIS at a global scale, especially in winter.

Influence of Snow on the Magnitude and Seasonal Variations of CIMODIS at a Global Scale
We summarize the temporal sequence of the average CIMODIS in snow-covered and snow-free areas for typical IGBP classes in the northern hemisphere at latitudes ≥ 30°N (Figure 4) to investigate the influence of snow on the magnitude and seasonal variation of CIMODIS. CIMODIS statistics for snow-free areas are presented in Table 2. For snow-covered areas, only the average CIMODIS in winter (from December to February) are presented in Figure 4 to ensure there are enough snow-covered pixels. The average CIMODIS of snow-covered areas for the EBF class are not presented in Figure 4 since very few pixels in the EBF class are covered by snow ( Figure 3). The average CIMODIS for EBF is calculated based on the pixels in the northern hemisphere at latitudes ≥ 0°N (instead of latitudes ≥ 30°N) since most EBF (>95%) are distributed in tropical and subtropical regions. Figure 4 shows that CIMODIS retrievals from the snow-covered pixels are considerably larger than the snow-free pixels. In addition, the biases of the CI retrieved by the main algorithm between the snow-covered and snow-free areas for open shrubs and herbaceous plants are larger than the biases for forests and closed shrubs since the open shrubs and herbaceous plants are more open than the forests and closed shrubs. A more open landcover type would have a stronger snow homogenization effect, which increases the bias of the CI between the snow-covered and snow-free areas. If the snow-covered pixels are not excluded from the statistics, a considerably larger seasonal variation will be observed for CIMODIS, even for the ENF class. However, the CIMODIS of the ENF class should be seasonally invariant since the structure of the ENF is relatively stable throughout the year [4]. Therefore, the CIMODIS retrieved from snow-covered pixels should be corrected or excluded from the statistics. It can be found from Figure 4 and Table 2 that the CIMODIS of different IGBP classes shows similar seasonal tendencies over the years. For deciduous vegetation classes, the clumping effect usually reaches a maximum in the summer and a minimum in the winter [39]. The clumping effect of evergreen forests, especially in EBF, is relatively stable throughout the year. In contrast, CIMODIS of the DBF class displays a large seasonal variation that varies from 0.67 in the summer to 0.81 in the winter. The magnitude of the seasonal variation for the DBF class is nearly 3 times that of  Table 2.
As shown in Figure 2, most CI MODIS values for snow-covered areas are either larger than 1.00 or have no value. Statistics indicate that in areas with latitudes ≥ 30 • N in the northern hemisphere, 94% of the pixels with a CI MODIS ≥ 1.00 are covered by snow. The snow on the ground is highly isotropic scattered, thus reducing the contrast between the sunlit and shadowed components of the vegetation canopy [58][59][60], resulting in a smaller NDHD. Therefore, a larger CI MODIS will be retrieved from the snow-covered pixels. The results further demonstrate that snow cover reduces the quality of the retrieved CI MODIS, resulting in either abnormally large CI values or no values. Therefore, the CI MODIS retrieved from snow-covered pixels should be corrected or excluded for various applications when studying the seasonal variation of CI MODIS at a global scale, especially in winter.

Influence of Snow on the Magnitude and Seasonal Variations of CI MODIS at a Global Scale
We summarize the temporal sequence of the average CI MODIS in snow-covered and snow-free areas for typical IGBP classes in the northern hemisphere at latitudes ≥ 30 • N ( Figure 4) to investigate the influence of snow on the magnitude and seasonal variation of CI MODIS . CI MODIS statistics for snow-free areas are presented in Table 2. For snow-covered areas, only the average CI MODIS in winter (from December to February) are presented in Figure 4 to ensure there are enough snow-covered pixels. The average CI MODIS of snow-covered areas for the EBF class are not presented in Figure 4 since very few pixels in the EBF class are covered by snow (Figure 3). The average CI MODIS for EBF is calculated based on the pixels in the northern hemisphere at latitudes ≥ 0 • N (instead of latitudes ≥ 30 • N) since most EBF (>95%) are distributed in tropical and subtropical regions. Figure 4 shows that CI MODIS retrievals from the snow-covered pixels are considerably larger than the snow-free pixels. In addition, the biases of the CI retrieved by the main algorithm between the snow-covered and snow-free areas for open shrubs and herbaceous plants are larger than the biases for forests and closed shrubs since the open shrubs and herbaceous plants are more open than the forests and closed shrubs. A more open landcover type would have a stronger snow homogenization effect, which increases the bias of the CI between the snow-covered and snow-free areas. If the snow-covered pixels are not excluded from the statistics, a considerably larger seasonal variation will be observed for CI MODIS , even for the ENF class. However, the CI MODIS of the ENF class should be seasonally invariant since the structure of the ENF is relatively stable throughout the year [4]. Therefore, the CI MODIS retrieved from snow-covered pixels should be corrected or excluded from the statistics. It can be found from Figure 4 and Table 2 that the CI MODIS of different IGBP classes shows similar seasonal tendencies over the years. For deciduous vegetation classes, the clumping effect usually reaches a maximum in the summer and a minimum in the winter [39]. The clumping effect of evergreen forests, especially in EBF, is relatively stable throughout the year. In contrast, CI MODIS of the DBF class displays a large Remote Sens. 2018, 10, 1194 7 of 15 seasonal variation that varies from 0.67 in the summer to 0.81 in the winter. The magnitude of the seasonal variation for the DBF class is nearly 3 times that of the EBF class. As shown in Table 2, forests are usually more clumped than shrubs and herbaceous plants. In addition, needleleaf forests tend to have smaller CI MODIS than broadleaf forests. The ENF class has the highest clumping effect, with an average CI of~0.56. In contrast, grassland (GL) is the least clumped class, with an average CI of~0.83.    and "Main-All" indicate CI retrievals by the main algorithm for snow-free, snow-covered, and both snow-free and snow-covered areas, respectively. Figure 5 shows a scatterplot that compares the CI retrieved by the main algorithm for the snow-free areas (CImain-snow-free) with the CI retrievals by the main algorithm (red dots, CImain-snow) or backup algorithm (blue dots, CIbackup-snow) for snow-covered areas. Figures 4 and 5 show that CImain-snow is considerably larger than CImain-snow-free with a bias of ~0.13 for needleleaf forests and ~0.32 for other vegetation types. CIbackup-snow shows better agreement with CImain-snow-free with a smaller RMSE and bias. The CIbackup-snow of the forests and herbaceous plants slightly underestimate CImain-snow-free, while the CIbackup-snow of the shrubs slightly overestimates CImain-snow-free. The magnitude of CIbackup-snow seems very stable because CI retrievals by the backup algorithm are calculated from the BRDF archetypes, which are the average BRDF shapes [39]. In summary, the backup algorithm can correct for the influence of snow to obtain a more reasonable CI for snow-covered areas. The backup algorithm can provide superior output in snow seasons probably because the CIs in the LUT of the backup algorithm are calculated based on the high-quality BRDF shapes of snow-free pixels. Thus, the CI calculated from the typical BRDF archetype in the LUT can represent the average CI of snow-free pixels for the corresponding AFX interval and should be closer to the real CIs than the CIs retrieved by the main algorithm in snow seasons.
The invariant CI retrieved by the backup algorithm in snow seasons is probably due to the relatively flat BRDF shapes of snow-covered pixels in snow seasons. We summarize the statistics of the AFX in the snow-covered areas on DOY 041, 2006 ( Figure 6). Figure 6 shows that the AFX values of snow-covered pixels are mainly distributed in the AFX interval of [0.87, 1.05]. This statistics further demonstrate that snow-covered pixels capture a narrow range of AFX variations around unity (i.e.,  and "Main-All" indicate CI retrievals by the main algorithm for snow-free, snow-covered, and both snow-free and snow-covered areas, respectively. Figure 5 shows a scatterplot that compares the CI retrieved by the main algorithm for the snow-free areas (CI main-snow-free ) with the CI retrievals by the main algorithm (red dots, CI main-snow ) or backup algorithm (blue dots, CI backup-snow ) for snow-covered areas. Figures 4 and 5 show that CI main-snow is considerably larger than CI main-snow-free with a bias of~0.13 for needleleaf forests and~0.32 for other vegetation types. CI backup-snow shows better agreement with CI main-snow-free with a smaller RMSE and bias. The CI backup-snow of the forests and herbaceous plants slightly underestimate CI main-snow-free , while the CI backup-snow of the shrubs slightly overestimates CI main-snow-free . The magnitude of CI backup-snow seems very stable because CI retrievals by the backup algorithm are calculated from the BRDF archetypes, which are the average BRDF shapes [39]. In summary, the backup algorithm can correct for the influence of snow to obtain a more reasonable CI for snow-covered areas. The backup algorithm can provide superior output in snow seasons probably because the CIs in the LUT of the backup algorithm are calculated based on the high-quality BRDF shapes of snow-free pixels. Thus, the CI calculated from the typical BRDF archetype in the LUT can represent the average CI of snow-free pixels for the corresponding AFX interval and should be closer to the real CIs than the CIs retrieved by the main algorithm in snow seasons.
The invariant CI retrieved by the backup algorithm in snow seasons is probably due to the relatively flat BRDF shapes of snow-covered pixels in snow seasons. We summarize the statistics of the AFX in the snow-covered areas on DOY 041, 2006 ( Figure 6). the AFX in the snow-covered areas on DOY 041, 2006 ( Figure 6). Figure 6 shows that the AFX values of snow-covered pixels are mainly distributed in the AFX interval of [0.87, 1.05]. This statistics further demonstrate that snow-covered pixels capture a narrow range of AFX variations around unity (i.e., [0.87, 1.05]). Thus, the CI values retrieved by the backup algorithm are almost invariant in the AFX interval of [0.87, 1.05] compared to the entire AFX interval of [0.5, 1.5] according to the CI-AFX LUT in Jiao et al. [39].

Discussion
In this study, we investigate the influence of snow on CIMODIS based on in situ CI measurements

Discussion
In this study, we investigate the influence of snow on CIMODIS based on in situ CI measurements and spaceborne CI product statistics. For in situ CI measurements, it is difficult to measure the CI on snow-covered days due to various environmental limitations, such as the safety of surveyors, the accessibility of the research area and the performance of instruments at low temperatures during

Discussion
In this study, we investigate the influence of snow on CI MODIS based on in situ CI measurements and spaceborne CI product statistics. For in situ CI measurements, it is difficult to measure the CI on snow-covered days due to various environmental limitations, such as the safety of surveyors, the accessibility of the research area and the performance of instruments at low temperatures during winter. Therefore, few CI measurements are available for snow-covered days in the current accumulated field-measured CI database. The field-measured CI values of snow-covered ENF in TP39 and TP74 [37] provide an opportunity to investigate the influence of snow on CI MODIS with real validation data. The results in Figure 1 indicate that CI retrievals using only the NDHD-CI equations usually present abnormally larger fluctuations than field-measured CIs. Similar results with considerably large CIs can also be found in the CI product derived by Wei and Fang [38] in the VALERI sites. In comparison, the CIs retrieved by the backup algorithm show good consistency with the field-measured CIs once such abnormal CI values occur in the main algorithm. Therefore, the results here show the efficiency of the backup algorithm in processing snow-covered pixels.
As indicated in the Figure 1, the CI MODIS retrieved by the main algorithm fluctuates too much with comparison to the in-situ measurements in snow-free seasons. It is still challenging to capture the seasonality of CI using spaceborne CI data since the seasonal variation of the CI may be contaminated by short-term fluctuations in the retrieved CI time series data [37]. Abnormally small or large values can also be found in the recent seasonal CI retrievals from MODIS BRDF products based on the CI-NDHD equations at validation sites [37][38][39]. In this study, we mainly focus on correcting the abnormally large CI values caused by the ephemeral snow in snow seasons. Further research is definitely needed to improve our understanding in the seasonal variation of CI, especially in rainy seasons.
For current available spaceborne CI products, CI are mainly retrieved based on the NDHD-CI equations developed by Chen et al. [4]. The cover of snow reduces the contrast between the sunlit and shadowed components of the vegetation canopy. Thus, the CI retrievals based on the NDHD-CI equations are greatly influenced by snow. The CI values of snow-covered pixels are usually marked as low quality and excluded from research when studying the seasonality of the CI [5,6,37,39], which may be appropriate depending on the research purpose. However, when the global CI product is used as the input parameter for ecological, hydrological and land surface models, the accuracy of the models will be greatly influenced if the CI values in the snow-covered areas are not corrected because more than 30 percent of the pixels in the northern hemisphere are covered by snow in winter. As shown in Figures 1, 4 and 5, the backup algorithm provides a good correction for the influence of snow to obtain more reasonable CI values. Therefore, in the latest version of our global CI product, the backup algorithm will be used to process the pixels marked as "snow" to further improve the quality of the CI product. After such a process, the CI product should be more suitable as the input parameter for ecological, hydrological and land surface models.

Conclusions
The CI and its seasonality are critical for global LAI estimating and ecological modeling. In this study, we investigate the influence of snow on the magnitude and seasonal variation of the CI retrieved from MODIS BRDF products based on the ground measurements of CI time series data for the three validation sites and global CI product statistics. The result highlights are as follows: (1) For snow-covered areas, the cover of snow will reduce the accuracy of CI retrievals by the main algorithm and result in abnormally larger CI values than the CI values for snow-free areas, with an average bias of~0.13 for needleleaf forests and~0.32 for other vegetation types. However, the backup algorithm can correct for the influence of snow to obtain more reasonable CI. Therefore, in the latest version of our global CI product, the backup algorithm will be further used to process the pixels of snow-covered areas. The CI product will have improved accuracy in snow-covered areas and is more suitable as the input parameters for ecological and meteorological models. (2) The time-series (e.g., 8-day) CI product shows the potential in investigating the trajectories of the clumping effect in snow seasons. However, extensive validation of the time-series CI product using field measurements is still a challenge, particularly in snow seasons. (3) The clumping effect of evergreen forests, especially for EBF, is relatively stable throughout the year. In contrast, the CI displays significant seasonal variation for deciduous vegetation types, particularly for DBF.
Snow cover information can not only be regarded as noise in CI determination but also as a clue to help retrieve a CI with higher accuracy. A more reasonable CI in snow-covered areas will be retrieved if the detailed information of snow (e.g., snow aging and ice re-crystallization processes) can be obtained from the corresponding product. In our future work, we will try to make full use of snow information to produce a global CI product that is more suitable as the input parameter for ecological and meteorological models.