## 2 Materials and Methods
### Retrieval Algorithm of the MODIS CI  
The CI retrieval algorithm using the hotspot-adjusted BRDF model [70] was exploited to calculate different versions of MODIS CIs. This algorithm framework includes the main algorithm and backup algorithm and constrains the retrieved CIs to the closed interval [0.33, 1.0] according to the theory of BRDF variability for vegetation canopy in the principal plane and, based on the relationship of the NDHD and CI. There are three main steps to retrieving the CI from MODIS multiangle reflectance data. First, the hotspot-adjusted BRDF model with two optimized hotspot parameters (C1 and C2) was used to fit the BRFs. Then, the reflectances in the direction of the hotspot (45\({}^{\circ}\) in backward scattering) and dark spot (45\({}^{\circ}\) in forward scattering) were reconstructed to calculate the NDHD. Third, the CI was retrieved in conjunction with IGBP data based on the CI-NDHD linear relationships established by the 4-scale model [24; 71; 72]. Each retrieved CI value was assigned a CI quality flag that is closely associated with the MODIS BRDF quality scheme. In particular, if the derived CI is outside the closed interval [0.33, 1.0] regardless of the high-quality flag from the MODIS BRDF product, a backup algorithm based on the AFX is adopted to process the so-called outliers.  
As the hotspot effect of surface anisotropic reflectance tends to be underestimated by the MODIS operational RTLSR model, a method has been developed to improve the estimation of the hotspot reflectance [73; 74; 75]. This method exploits an exponential approximation of the hotspot kernel to correct the hotspot effect of the volumetric scattering component for the RTLSR model. Two prior hotspot parameter values (i.e., C1 = 0.7, C2 = 3.2) were used in the improved model (i.e., the hotspot-adjusted version of the RTLSR model (RTCLSR)) to reconstruct the reflectance of hotspots (\(\rho_{HS}\)) and dark spots (\(\rho_{DS}\)). Then, the NDHD can be obtained by the following:  
\[\text{NDHD}=\frac{\rho_{HS}-\rho_{DS}}{\rho_{HS}+\rho_{DS}} \tag{2}\]  
According to the linear CI-NDHD equations in the red band [24; 71], the CI can be obtained in association with _a priori_ classification information. For coniferous forests,  
\[\text{CI}=-0.47\times\text{NDHD}+0.80 \tag{3}\]For the other vegetation classes,  
\[\text{CI}=-1.23\times\text{NDHD}+1.34 \tag{4}\]  
For the backup algorithm, the AFX was used to find the most similar BRDF archetype of the same land cover type and provide a most likely estimation of the CI for pixels outside the close interval [0.33, 1.0] in terms of the main algorithm. The AFX is calculated using the BRDF model parameters as follows:  
\[\text{AFX}=\frac{WSA(\lambda)}{f_{iso}(\lambda)}=1+\frac{f_{rol}(\lambda)}{f_{ iso}(\lambda)}\times H_{vol}+\frac{f_{geo}(\lambda)}{f_{iso}(\lambda)}\times H_{ geo} \tag{5}\]  
where \(WSA(\lambda)\) is the white sky albedo and \(H_{vol}\) and \(H_{geo}\) are bi-hemispherical integral values of volumetric and GO scattering kernels, respectively. The corresponding CI value of the backup algorithm can be found through a CI-AFX look-up table (LUT) established by using the CI-NDHD relationship based on the average BRDF shape (i.e., BRDF archetype) determined by the AFX.  
To analyze the C5 and C6 CI products on a monthly scale, temporal compositing is utilized to select the best retrievals and generate monthly CIs from 8-day (C5) and daily (C6) retrievals. Compositing mainly consists of two steps. The first step is selecting retrievals according to quality flags: the main algorithm retrievals with the best quality MCD43A1 data have the highest priority, followed by the main algorithm retrievals with good-quality MCD43A1 data, and if no main algorithm retrievals are available, backup algorithm retrievals are selected. Then, the monthly CI is calculated as the average of the selected retrievals, and corresponding quality flags are used to identify monthly QA data.